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A  COMPARISON  OF  SEVERAL  COASTAL  OCEAN  MODELS 


1.0  INTRODUCTION 

This  report  presents  results  from  a  comparison  of  several  ocean  models  that  are  being  used  for 
coastal  simulation  and  prediction.  The  models  are  tested  to  determine  their  ability  to  simulate  basic 
physical  processes  of  importance  in  the  coastal  ocean.  The  purpose  of  the  comparison  is  to  look  at 
the  pros  and  cons  of  the  models,  to  identify  particular  problems,  and  to  provide  a  baseline  for 
testing  future  models  and  making  improvements.  The  overall  goal  is  to  try  to  determine  the  best 
models  and/or  parameterizations  for  either  general  or  particular  applications. 

Physical  processes  important  in  the  coastal  ocean  include  tides,  wind  setup  and  storm  surge, 
wind-driven  circulation,  river  outflows  and  fresh  water  runoff,  surface  and  internal  waves, 
advection,  mixing,  coastal-trapped  waves,  upwelling  and  downwelling,  and  frontal  dynamics.  With 
the  exception  of  short-spatial-scale,  non-hydrostatic  effects,  e.g.,  such  as  occur  with  short- 
wavelength,  high-amplitude  internal  waves,  these  processes  can  generally  be  adequately  simulated 
with  hydrostatic,  three-dimensional  (3-D)  models  of  the  type  to  be  compared  in  this  study. 

We  note  that  sediment  and  biological  processes  and  their  effect  on  the  water’s  optical  properties 
are  also  very  important  in  the  coastal  environment.  However,  investigation  of  these  processes, 
which  requires  that  sediment,  biological,  and  optical  submodels  be  linked  to  the  ocean  circulation 
models,  is  not  addressed  here. 

The  model  comparison  was  restricted  to  ocean  models  that  include  a  free  surface  and  predict 
the  3-D  fields  of  ocean  currents,  temperature,  and  salinity.  A  number  of  widely  used  ocean  models 
employ  a  rigid  lid  to  filter  out  surface  gravity  waves  and  avoid  the  timestep  restriction  imposed  by 
surface  waves.  Rigid-lid  models  can  be  used  to  investigate  many  important  coastal  processes. 
However,  a  free  surface  is  needed  for  the  prediction  of  tides  and  storm  surge.  Hence,  a  model  that 
is  to  be  used  to  simulate  the  complete  range  of  coastal  processes  should  have  a  free  surface. 

The  models  compared  here  include  the  Princeton  Ocean  Model  (POM),  which  was  originally 
developed  by  Alan  Blumberg  and  George  Mellor  (Blumberg  and  Mellor  1987);  the  Estuarine  and 
Coastal  Ocean  Model,  Semi-Implicit  version  (ECOM-si)  of  Alan  Blumberg,  which  is  based  on 
the  original  Princeton  model,  but  has  some  significant  modifications  (Blumberg  1992);  the  Sigma 
z-level  model  (SZM),  developed  in  house  at  the  Naval  Research  Laboratory  (NRL)  (Martin  1998); 
and  the  S-Coordinate  Rutgers  University  Model  (SCRUM),  which  was  developed  by  Tony  Song  and 
Dale  Haidvogel  (Song  and  Haidvogel  1994). 

POM  is  probably  the  most  widely  used  of  all  baroclinic  coastal  ocean  models  and  has  been 
applied  to  a  wide  range  of  coastal  problems.  A  list  of  publications  in  which  POM  has  been  used 
is  available  from  the  POM  ftp  site  at  Princeton  University  (ftp.gfdl.gov)  and  has  over  100  listings. 
POM  was  originally  developed  by  Alan  Blumberg  and  George  Mellor  in  the  early  1980s.  Since  that 
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time,  there  have  been  contributions  to  the  program  by  others,  notably  Leo  Oey,  Boris  Galperin,  and 
Lakshmi  Kantha.  The  major  features  of  POM  are  (a)  the  use  of  an  Arakawa  “C”  grid,  (b)  second- 
order,  centered  spatial  finite  differencing,  (c)  leapfrog  time  differencing  with  an  Asselin  temporal 
filter  to  eliminate  time  splitting,  (d)  an  explicit  treatment  of  the  surface  waves  using  a  smaller 
timestep  than  that  used  for  the  internal  mode,  (e)  a  sigma  (bottom  following)  vertical  coordinate, 
(f)  an  implicit  treatment  of  vertical  mixing,  and  (g)  the  use  of  the  Mellor-Yamada  Level  2.5 
(MYL2.5)  turbulence  scheme  to  calculate  vertical  mixing. 

ECOM-si  is  based  on  POM,  but  contains  some  significant  modifications  that  were  implemented 
by  Alan  Blumberg,  Vincenzo  Casulli,  and  Ralph  Cheng.  The  main  differences  from  POM  are  (a) 
the  use  of  a  two-time-level  temporal  scheme  rather  than  leapfrog,  (b)  the  use  of  an  implicit  rather 
than  an  explicit  scheme  for  the  free  surface,  and  (c)  the  addition  of  a  wetting/drying  capability. 

SZM  was  developed  at  NRL  to  provide  a  model  than  can  use  both  sigma  and  z-level  (i.e.,  fixed 
depth)  vertical  coordinates.  The  model  uses  sigma  coordinates  down  to  a  user-specified  depth  and 
z-levels  below.  This  allows  some  flexibility  in  setting  up  the  vertical  grid  and  allows  comparisons 
to  be  made  between  sigma  and  z-level  coordinates.  In  other  respects,  SZM  is  similar  to  POM, 
except  that  it  has  an  implicit  treatment  of  the  free  surface  and  uses  the  simpler  Mellor-Yamada 
Level  2  (MYL2)  turbulence  scheme. 

SCRUM  is  one  of  a  variety  of  ocean  models  that  have  been  developed  by  Dale  Haidvogel’s 
ocean  modeling  group  at  Rutgers  University.  Dale  Haidvogel’s  group  is  known  for  trying  a  wide 
range  of  ocean  modeling  techniques,  including  spectral  and  finite  element  (FE)  methods.  The  best 
known  of  their  models  is  the  Semi-spectral  Primitive  Equation  Model  (SPEM)  (Haidvogel  et  al. 
1991).  SPEM  has  been  used  in  a  number  of  modeling  studies  that  have  appeared  in  the  oceanographic 
literature;  however,  SPEM  is  limited  for  coastal  applications  because  of  its  rigid  lid. 

SCRUM  (Version  2.1)  uses  finite  differences  (FD)  in  the  horizontal,  but  uses  FE  in  the  ver¬ 
tical.  Other  major  differences  from  the  other  models  being  compared  here  are  (a)  the  use  of  a 
third-order  Adams-Bashforth  temporal  scheme  for  most  the  baroclinic  terms,  (b)  the  use  of  the 
Crank-Nicolson  semi-implicit  scheme  for  vertical  advection  and  mixing,  and  (c)  the  use  of  a 
generalized  sigma  coordinate  system  in  the  vertical,  called  an  “S”  coordinate  by  Song  and  Haidvogel, 
that  allows  changing  the  relative  spacing  of  the  vertical  layers  as  the  depth  of  the  water  changes. 
The  “S”  coordinate  allows  greater  flexibility  in  setting  vertical  resolution  than  a  standard  sigma 
coordinate,  e.g.,  it  allows  one  to  maintain  a  minimum  resolution  in  the  surface  and/or  bottom 
boundary  layers  as  the  water  depth  increases.  SCRUM  uses  a  split-explicit  scheme  for  the  free 
surface  and  offers  a  choice  of  turbulence  schemes,  MYL2,  Price,  and  Pacanowski  and  Philander  (PP) 
to  calculate  the  vertical  mixing  coefficients. 

A  problem  that  was  encountered  in  testing  SCRUM  was  that  during  the  period  over  which  the 
model  comparisons  were  conducted,  SCRUM  was  still  undergoing  some  development.  The  version 
of  SCRUM  that  was  started  with,  Version  2.1,  was  discovered  to  have  some  problems,  and  testing 
was  begun  on  the  other  three  models  while  waiting  for  a  more  fully  developed  version  of  SCRUM  to 
become  available.  About  this  time,  Tony  Song  (who  did  most  of  the  original  development  of 
SCRUM)  left  Rutgers,  and  the  development  of  SCRUM  at  Rutgers  was  continued  by  Hernan 
Arango.  SCRUM  Version  3.0  was  released  by  Rutgers  (Arango,  pers.  comm.)  in  August  1996.  This 
version  of  SCRUM  differed  significantly  from  the  previous  version,  a  major  change  being  that  the 
FE  scheme  used  in  the  vertical  had  been  replaced  by  a  FD  scheme  similar  to  that  used  by  the  other 
models.  Because  of  the  goal  of  finishing  up  the  initial  phase  of  the  model  comparison  study  at  this 
time,  it  was  deemed  too  late  to  begin  testing  a  new  version  of  SCRUM.  Hence,  the  discussion  of 


A  Comparison  of  Several  Coastal  Ocean  Models 


3 


SCRUM  in  this  report  refers  to  Version  2.1  that  was  originally  received.  Since  a  fully  sorted 
version  of  this  model  was  never  obtained,  most  of  the  tests  that  were  conducted  do  not  include 
results  from  SCRUM. 

Besides  the  models  being  tested  here,  there  are,  of  course,  a  significant  number  of  other  coastal 
models  in  current  use.  Practical  considerations  limited  the  number  of  models  that  could  be  included 
in  this  study,  but  some  particular  models  of  interest  will  be  briefly  mentioned. 

Curvilinear  Hydrodynamics  in  3  Dimensions  (CH3D)  is  a  coastal,  baroclinic  model  that  is 
being  used  by  the  Coastal  Engineering  Research  Center  (CERC)  of  the  U.  S.  Army  Corps  of 
Engineers  (Johnson  et  al.  1991).  This  model  was  initially  developed  for  CERC  by  Peter  Sheng 
of  the  University  of  Florida,  who  has  been  active  in  coastal  and  lake  modeling  since  the  1970s. 
CERC  is  a  FD  model  and  is  similar  in  formulation  to  POM.  In  recent  years,  CH3D  has  been 
maintained  and  developed  by  CERC. 

The  Tidal,  Residual,  Intertidal  Mudflat,  Three-Dimensional  Model  (TRIM-3D)  is  a  FD  coastal 
and  estuary  model  developed  by  Vince  Casulli  and  Ralph  Cheng  (Casulli  and  Cheng  1993).  A 
unique  feature  of  TRIM-3D  is  that  grid  cells  over  land  are  not  included  in  the  model’s  storage 
arrays  or  calculations.  This  allows  a  very  large  increase  in  efficiency  in  simulating  coastal  and 
estuary  systems  in  which  a  large  fraction  of  the  volume  encompassed  by  the  domain  being  inves¬ 
tigated  is  land.  TRIM-3D  has  been  applied  to  estuary  systems  in  which  the  ratio  of  sea  points  to 
land  points  within  the  3-D  volume  encompassing  the  estuary  system  is  less  than  2%  (Casulli  and 
Cheng  1993).  NRL  has  recently  developed  a  version  of  POM  in  which  land  areas  have,  in  a  similar 
fashion,  been  eliminated  from  the  array  storage  and  model  calculations  (Ko,  pers.  comm.). 

David  Dietrich  of  Mississippi  State  University  has  developed  a  free  surface  version  of  his 
DieCAST  model  (Dietrich  and  Mehra  1998).  The  original  version  of  DieCAST,  which  uses  a  rigid 
lid,  is  able  to  form  and  maintain  mesoscale  circulation  features  and  fronts  with  relatively  low 
horizontal  resolution  (Dietrich  et  al.  1993;  Dietrich  and  Ko  1994). 

Dan  Lynch  of  Dartmouth  University,  one  of  the  pioneers  in  the  development  and  application 
of  barotropic  FE  ocean  models,  has  developed  a  baroclinic  FE  model  called  QUODDY3  (Lynch  and 
Werner  1991).  This  model  is  currently  being  applied  by  Dan  Lynch  and  colleagues  to  the  Maine 
Coastal  Current  (Naimie  et  al.  1994).  FE  models  are  more  complex  to  program  than  FD  models, 
but  their  very  flexible  horizontal  grids  allow  large  changes  in  spatial  resolution  over  short  distances, 
which  provides  significant  advantages  in  modeling  coastal  regions  with  complex  coastlines  and 
bathymetry.  (NRL  has  recently  acquired  and  begun  working  with  QUODDY3.) 

The  models  being  compared  in  this  study,  POM,  ECOM-si,  SZM,  and  SCRUM,  although 
similar  in  many  ways,  have  some  significant  differences:  (a)  explicit  versus  implicit  treatment  of 
the  free  surface,  (b)  sigma  versus  z-level  vertical  coordinates,  (c)  different  temporal  differencing 
schemes,  and  (d)  different  vertical  mixing  submodels.  This  study  will  consider  how  these  differences 
affect  model  performance. 

The  model  tests  conducted  here  consist  of  tests  of  basic  physical  processes,  including  advection, 
mixing,  propagation  of  free  and  coastal-trapped  waves,  and  formation  of  upwelling  and  down- 
welling  fronts.  There  are  several  reasons  for  conducting  tests  of  basic  physical  processes,  rather 
than  tests  of  more  complex  coastal  situations  and  comparisons  with  observations: 

•  Coastal  models  are  sometimes  applied  without  good  knowledge  of  how  well  they  simulate 
basic  physical  processes. 
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•  Knowledge  of  model  performance  in  simulating  basic  processes  can  aid  interpretation  of 
results  in  realistic  situations  and  help  in  identifying  the  cause  of  model-related  problems. 

•  Model  comparisons  with  real  data  are  not  always  conclusive  because  of  uncertainties  in  initial 
conditions,  boundary  conditions,  forcing,  and  validation  data. 

•  Particular  simulations  with  real  data  may  not  provide  a  good  test  of  all  the  important  physical 
processes. 

•  Tests  of  basic  processes  can  uncover  problems  with  the  models  that  may  not  be  evident  in 
particular  simulations  with  real  data. 

It  is  acknowledged,  however,  that  comparison  with  observations  is  the  final  arbiter  of  model 
skill.  There  are  a  number  of  studies  being  conducted  with  coastal  ocean  models  at  NRL  and  at  other 
institutions  where  the  primary  focus  is  the  comparison  of  model  results  with  observations.  We 
consider  the  model  tests  conducted  here  to  be  complimentary  to  these  other  studies. 

The  sections  that  follow  include  (2.0)  a  description  of  the  models,  (3.0)-(8.0)  a  discussion  of 
the  results  of  the  various  model  comparison  test  cases,  and  (9.0)  a  summary. 


2.0  DESCRIPTION  OF  THE  MODELS 

The  description  of  the  models  provided  here  is  not  complete  in  all  details.  The  reader  is 
referred  to  the  references  provided  for  the  models  for  a  more  complete  description.  The  emphasis 
here  is  on  how  the  model  formulations  differ  from  each  other.  Hence,  rather  than  sequentially 
present  a  complete  description  of  each  model,  the  different  aspects  of  the  models,  i.e.,  the  horizontal 
and  vertical  coordinate  systems  used,  the  spatial  and  temporal  differencing,  the  treatment  of  the  free 
surface,  and  the  parameterization  of  horizontal  and  vertical  mixing,  are  discussed  in  turn  for  all 
the  models.  There  is  some  additional  discussion  of  the  physics  and  numerics  of  the  models  in  the 
sections  describing  the  model  tests.  A  listing  of  some  of  the  main  features  of  the  models  is 
presented  in  Table  1. 


2.1  Basic  Equations 

All  the  models  have  a  free  surface  and  employ  the  hydrostatic,  Boussinesq,  and  incompressible 
approximations.  The  equations  that  the  models  solve,  written  in  Cartesian  coordinates,  are 


P  =  ~?8 


du  dv  dw 
dx  +  dy  +  dz 


(1) 

(2) 

(3) 


(4) 
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Table  1  —  Comparison  of  Some  of  the  Features  of  Models  Being  Tested  in 

Comparison  Study 


POM 

ECOM-si 

SZM 

SCRUM 

Horizontal  Grid 

C-Grid 

C-Grid 

C-Grid 

C-Grid 

Curvilinear 

Curvilinear 

Cartesian 

Curvilinear 

Vertical  Grid 

Sigma 

Sigma 

Sigma/z-  level 

S-Coord 

Barotropic  Mode 

Split-Explict 

Leapfrog 

Implict 

Two-Time  Level 

Semi-Explict 

Leapfrog 

Split-Explicit 

Leapfrog 

Horizontal  FD 

Second-Order 

Second-Order 

Second-Order 

Second-Order 

Vertical  FD 

Second-Order 

Second-Order 

Second-Order 

Linear  FE 

Horizontal  Mixing 

Laplacian 

Smagorinsky 

Laplacian 

Smagorinsky 

Laplacian 

Grid  Cell-RE 

Laplacian  or 
Biharmonic 

Vertical  Turbulence 

MYL2.5 

MYL2.5 

MYL2 

MYL2 

Price  PP 

Bottom  Friction 

Quadratic 

Quadratic 

Quadratic 

Quadratic 

or  Linear 

or  Linear 

Wetting/Drying 

No 

Yes 

No 

No 

Temporal  Scheme  for  Barotropic  Equations 

Surface  Gradient 

Centered 

Fully  Implicit 

Trapezoidal 

Centered 

Transport  Gradient 

Centered 

Fully  Implicit 

Trapezoidal 

Centered 

Temporal  Scheme  for  Baroclinic  Equations 

Pressured  Gradient 

Centered 

Forward 

Centered 

Adams-Bashforth 

Coriolis  Term 

Centered 

Forward 

Centered 

Adams-Bashforth 

Horizontal  Advection 

Centered 

Forward 

Centered 

Adams-Bashforth 

Horizontal  Mixing 

Forward 

Forward 

Forward 

Adams-Bashforth 

Vertical  Advection 

Centered 

Forward 

Centered 

Crank-Nicolson 

Vertical  Mixing 

Fully  Implicit 

Fully  Implicit 

Fully  Implicit 

Crank-Nicolson 

dt_  d((t  +  ff)u)  d(tt  +  ff)v) 

dt  dx  dy  ^ 

f  -  -v '  <vr> + v*  <w> + I  <*«!> + Q'Tz  •  (6> 

§  =  -V  •  (VS)  +  V*  (A„V„S)  +  |  (K„£)  ,  (7) 

p  =  p(r,  s,  z ) ,  (8) 

where  t  is  the  time,  x,  y,  and  z  are  the  three  coordinate  directions,  u,  v,  and  w  are  the  three 
components  of  the  vector  velocity  v,  T  is  the  potential  temperature,  5  is  the  salinity,  Vh  is  the 
horizontal  gradient  operator,  /  is  the  Coriolis  parameter,  p  is  the  pressure,  p  is  the  water  density. 
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p0  is  a  reference  density,  g  is  the  acceleration  of  gravity,  Fu  and  Fv  are  horizontal  mixing  terms 
for  momentum,  AH  is  the  horizontal  mixing  coefficient  for  scalar  fields  (temperature  and  salinity), 
Km  and  Kh  are  vertical  eddy  coefficients  for  the  momentum  and  scalar  fields,  respectively,  Qr  is 
the  solar  radiation,  y  is  a  function  describing  the  solar  extinction,  £  is  the  elevation  of  the  free 
surface  above  the  undisturbed  value  at  z  =  0,  H  is  the  bottom  depth,  and  u  and  V  are  the 
depth-averaged  horizontal  velocities. 

With  the  hydrostatic  assumption,  vertical  accelerations  are  assumed  to  be  small  and  the  vertical 
momentum  equation  (3)  is  taken  to  be  a  balance  between  the  vertical  pressure  gradient  and  the 
gravitational  force.  Vertical  velocity  is  computed  from  the  divergence  of  the  horizontal  flow  field 
using  the  continuity  equation  (4). 

With  the  Boussinesq  approximation,  density  variations  are  only  taken  into  account  in  calculating 
the  horizontal  pressure  gradient  terms  and  in  calculating  vertical  stability  for  the  parameterization 
of  vertical  mixing. 

All  the  models  use  some  version  of  the  nonlinear  equation  of  state  for  seawater  to  calculate 
density  (8)  from  T  and  S.  ECOM-si,  SZM,  and  SCRUM  calculate  potential  density,  i.e.,  the  density 
at  atmospheric  pressure.  POM  calculates  the  in  situ  density,  which  includes  the  effect  of  pressure, 
to  provide  more  accurate  density  gradients  in  deep  water.  Note  that  the  in  situ  density  must  be 
corrected  for  the  effect  of  pressure  when  calculating  vertical  stability  (which  POM  does). 


2.2  Surface  and  Bottom  Boundary  Conditions 

Equations  (l)-(8)  are  subject  to  boundary  conditions  in  the  form  of  fluxes  and  stresses  at  the 
ocean’s  surface  and  bottom.  The  boundary  conditions  at  the  surface,  which  are  due  to  fluxes 
between  the  ocean  and  the  atmosphere,  are 


du  xx 

KuTz-y0' 

(9) 

dv  xy 

KmTz  =  T„’ 

(10) 

dT  Qb+Qe+  Qs 

Kfi  <)Z  poCp 

(11) 

K«%^S\Z,0(EV-Pr), 

(12) 

where  t*  and  xy  are  the  x  and  y  components  of  the  surface  wind  stress,  Qb,  Qe,  and  Qs  are  the  net 
long-wave,  latent,  and  sensible  surface  heat  fluxes,  Ev  and  Pr  are  the  surface  evaporation  and 
precipitation  rates,  and  cp  is  the  specific  heat  of  seawater.  At  the  ocean  bottom,  the  boundary 
conditions  are 


du  _  Tb 


KMx~  = 


(13) 
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Km 


dv  x 


&  Po  ’ 
dT 

K»ir-0- 

dS 

KhTz=0, 


(14) 


(15) 

(16) 


where  x*  and  x£  are  the  x  and  y  components  of  the  bottom  stress.  The  bottom  stress  is  parameterized 
either  by  a  linear  drag  law 


Tl  =  -P ocbxu  * 

(17) 

Xl  =  -P ocbxv  > 

(18) 

where  c ^  is  a  linear  drag  coefficient  with  units  of  velocity,  or  a  quadratic  law 


\  =-Pocfe«M» 


(19) 


xl  =  -PocfcvM  »  (2°) 

where  cb  is  a  dimensionless  coefficient. 

The  models  generally  use  a  quadratic  bottom  drag  law  with  the  drag  coefficient  computed  in 
terms  of  the  bottom-layer  thickness  A zb  and  the  bottom  roughness  z0  as 


cb  =  max 


2 

K 

log2 

v2z°) 

cb„ 


(21) 


where  k  =  0.4  is  von  Karman’s  constant  and  cb  is  a  minimum  value  for  cb.  This  expression  for 
cb  is  derived  by  combining  (19-20)  with  the  equation  for  a  logarithmic  boundary  layer  velocity 
profile 


(22) 


where  z  is  the  distance  above  the  bottom. 
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2.3  Horizontal  Coordinate 

POM,  ECOM-si,  and  SCRUM  use  general  orthogonal,  curvilinear,  horizontal  coordinates.  SZM 
is  written  to  use  only  Cartesian  coordinates  and  requires  the  use  of  a  constant  horizontal  grid 
spacing  in  x  and  y. 

Curvilinear  grids  allow  some  flexibility  in  setting  up  the  horizontal  grid  in  terms  of  varying  the 
horizontal  resolution  in  different  parts  of  the  model  domain  and  positioning  the  lateral  boundaries, 
or  alternatively,  conforming  to  a  particular  map  projection.  There  are,  however,  some  limitations 
on  the  degree  to  which  an  orthogonal,  curvilinear  grid  can  be  curved  or  the  rate  at  which  the  grid 
resolution  can  be  changed  if  truncation  errors  associated  with  the  spatial  changes  in  the  grid  are  to 
be  kept  small  (Crowder  and  Dalton  1971;  Rivas  1972). 

The  implementation  of  orthogonal  curvilinear  coordinates  in  the  models  is  primarily  a  matter 
of  accounting  for  the  changing  size  of  the  grid  cells  when  computing  fluxes  between  the  cells  and 
when  computing  the  changes  within  a  cell.  Hence,  the  variables  for  the  horizontal  size  of  the  grid 
cells,  Ax  and  Ay,  must  be  stored  as  two-dimensional  (2-D)  arrays  rather  than  as  constants.  The  only 
correction  to  the  models’  equations  for  the  curvature  of  the  horizontal  grid  is  a  correction  to  the 
advection  term  to  account  for  conversion  between  u  and  v  momentum  for  advection  along  curving 
grid  coordinates.  Note  that  the  horizontal  momentum  diffusion  terms  should  also,  ideally,  have  a 
curvature  correction.  Since  the  models  do  not  have  such  a  correction,  they  implicitly  assume  that 
the  diffusion  error  due  to  any  curvature  of  the  grid  is  small  and  can  be  neglected.  None  of  the 
models  account  for  changes  in  the  grid-cell  spacing  when  interpolating  variables  between 
the  center  and  the  interfaces  of  the  grid  cells,  i.e.,  all  the  models  use  simple  averages  for  these 
interpolations. 

All  of  the  tests  in  this  report  were  conducted  with  Cartesian  horizontal  coordinates  with  Ax  and 
Ay  both  constant.  Hence,  error  that  might  arise  from  the  particular  implementation  of  curvilinear 
horizontal  coordinates  in  the  individual  models  was  not  addressed. 


2.4  Vertical  Coordinate 

POM  and  ECOM-si  use  a  sigma  (a)  coordinate  in  the  vertical  in  which  the  depths  of  the  model 
layers  are  a  fixed  fraction  of  the  total  depth  of  the  water  column  from  the  free  surface  to  the  bottom. 
The  sigma  vertical  coordinate  is  expressed  as 


a  = 


ML 

(?+")’ 


(23) 


so  that  0  varies  from  0  at  the  surface  to  -1  at  the  bottom.  With  sigma  coordinates,  the  layers  near 
the  bottom  follow  the  contours  of  the  bottom. 

SZM  uses  sigma  coordinates  for  a  user-specified  number  of  upper  layers  and  uses  z-levels  for 
the  layers  below.  SZM  can  be  run  with  all  sigma  layers,  with  sigma  layers  in  shallow  water  and 
z-levels  in  deep  water,  or  with  all  z-levels  (except  for  the  upper  layer,  which,  because  of  the  free 
surface,  must  be  a  sigma  layer).  All  the  models  allow  stretching  of  the  vertical  grid,  i.e.,  changes 
in  Ao  or  A z  with  depth,  to  provide  different  vertical  resolution  over  different  parts  of  the  water 
column. 
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SCRUM  uses  what  Song  and  Haidvogel  (1994)  call  an  “S”  vertical  coordinate.  This  vertical 
coordinate  allows  changing  the  fractional  size  of  the  individual  sigma  layers  (Act  =  A z/(£  +  H))  in 
different  parts  of  the  domain  as  the  depth  of  the  bottom  changes,  rather  than  keeping  Act  for  a 
particular  layer  constant  everywhere,  as  is  the  case  with  the  usual  implementation  of  sigma  layers 
used  in  POM,  ECOM-si,  and  SZM.  For  example,  a  minimum  thickness  can  be  specified  for 
layers  near  the  surface  and  bottom  to  avoid  the  normal  vertical  spreading  of  the  sigma  layers  as  the 
depth  of  the  water  increases. 

Similar  to  orthogonal,  curvilinear  coordinates,  the  implementation  of  sigma  or  “S”  coordinates 
in  the  models  is  primarily  a  matter  of  accounting  for  the  changing  vertical  thickness  of  the  layers  in 
calculating  fluxes  between  adjacent  grid  cells  and  in  calculating  changes  within  the  grid  cells.  Note 
that  with  sigma  coordinates,  the  changes  in  the  thickness  of  the  grid  cells  occurs  not  only  horizontally 
within  a  layer,  but  also  from  timestep  to  timestep  because  of  the  changing  surface  elevation.  To 
avoid  having  to  use  a  large  number  of  3-D  arrays  to  store  the  thickness  of  the  grid  cells  at  different 
locations  and  time  levels,  the  models  express  the  grid  thickness  A z  as  the  product  of  the  fractional 
thickness  of  the  sigma  layer  A  a  times  the  total  depth,  i.e.,  as 

Az  =  Ao  (£  +  H)  .  (24) 

The  only  correction  to  the  equations  used  in  the  models  for  the  changing  depth  of  the  sigma  layers 
is  a  correction  for  the  horizontal  pressure  gradient.  The  horizontal  pressure  gradient  in  sigma 
coordinates  contains  an  extra  term  to  calculate  and  remove  the  vertical  change  in  pressure  between 
neighboring  points  within  a  sigma  layer  so  that  the  net  pressure  change  that  is  computed  will  be 
approximately  along  a  level  surface  (Blumberg  and  Mellor  1987). 

A  problem  with  this  calculation  of  the  pressure  gradient  with  sigma  coordinates  is  that  the 
vertical  component  of  the  pressure  change  along  a  sloping  sigma  surface  is  frequently  much  larger 
than  the  horizontal  component  (Haney  1991).  In  this  case,  the  desired  horizontal  component  is 
calculated  as  the  small  difference  between  two  large  terms  and  is  subject  to  significant  truncation 
error.  An  expedient  that  is  employed  in  the  models  to  reduce  this  error  is  to  subtract  the  horizontally 
averaged  density  profile  from  the  3-D  density  field  when  calculating  the  horizontal  pressure  gradient 
so  that  the  main  component  of  the  vertical  change  in  pressure  is  removed  from  the  calculation. 

Strictly  speaking,  for  a  full  transformation  of  the  equations  to  sigma  coordinates,  the  horizontal 
diffusion  terms  would  also  be  corrected  for  the  transformation  so  that  they  would  still  repre¬ 
sent  diffusion  along  level  surfaces.  However,  all  the  models  use  the  approximation  discussed  by 
Mellor  and  Blumberg  (1985)  who  argued  that  diffusion  along  the  sigma  surfaces  rather  than  along 
level  surfaces  was,  in  general,  more  appropriate  for  sigma  coordinate  models,  particularly  for 
proper  simulation  of  the  bottom  boundary  layer  with  sigma  coordinates.  A  practical  reason  for  having 
“horizontal”  diffusion  occur  along  the  sigma  surfaces  is  that  the  rather  messy  form  of  the  transformed 
horizontal  diffusion  terms  is  avoided. 


2.5  Spatial  Finite  Differences 

All  the  models  use  standard,  second-order,  centered  spatial  interpolations  and  FDs  except  that 
SCRUM  uses  linear  FEs  in  the  vertical.  With  second-order,  centered  interpolations  and  differences, 
the  value  of  a  variable  <J>  at  a  location  x  is  evaluated  as  the  average  of  the  values  defined  on  either 
side,  i.e., 
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♦lx 


1 

2 


<t>  1. 

x  +  -Ax 
2 


+  <t> 


and  the  gradient  of  (j)  at  x  is  calculated  as 


dx'X 


Ax 


4>  iA 

x  +  —Ax 
2 


(25) 


(26) 


where  Ax  is  the  grid  spacing. 


2.6  Temporal  Scheme  for  Baroclinic  Equations 

The  temporal  schemes  that  are  used  for  the  3-D  baroclinic  equations  of  the  models  are  discussed 
in  this  section.  The  temporal  schemes  are  illustrated  with  just  the  u  momentum  equation,  since  the 
treatment  of  the  other  model  variables  is  similar.  In  the  following  discussion,  ( n )  will  denote  model 
values  at  the  current  time  level  (i.e.,  values  calculated  on  the  previous  timestep),  (n  +  1)  will  denote 
the  newly  calculated  values,  and  («  -  1)  and  (n  -  2)  will  denote  values  at  previous  time  levels.  Note 
that,  for  simplicity,  the  temporal  differencing  equations  that  are  presented  here  do  not  include  the 
temporal  changes  in  the  vertical  thickness  of  the  layers  on  the  models’  sigma-coordinate  grids  (and 
on  SCRUM’s  S-coordinate  grid)  caused  by  changes  in  the  surface  elevation  (though  these  changes 
are  accounted  for  by  the  models). 

POM  and  SZM  use  a  standard,  leapfrog  time-stepping  scheme.  With  this  scheme,  most  of  the 
terms,  i.e.,  the  advection,  baroclinic  pressure  gradient,  and  Coriolis  terms,  are  centered  in  time  at 
(«).  The  horizontal  diffusion  terms  are  calculated  at  the  previous  (n  -  1)  time  level  (since  diffusion 
terms  evaluated  at  the  central  time  level  of  a  leapfrog  scheme  tend  to  excite  time  splitting  and  cause 
numerical  instability),  and  the  vertical  diffusion  terms  are  treated  implicitly  so  as  to  avoid  the 
timestep  restriction  for  explicit  vertical  diffusion  (the  high  rates  of  vertical  diffusion  sometimes 
calculated  by  the  turbulence  submodels  would  require  a  very  small  timestep  for  stability  if  the 
vertical  diffusion  were  explicit).  Hence,  the  temporal  differencing  of  the  u  momentum  equation  can 
be  represented  as 


m("  +  1)_m("-1)  , ,  , .  i  dDW  /„  n 

ar - v.(v„)<",+/v<”)-i^r+F<''-,)+ 


K, 


M ' 


dz 


r 


(27) 


POM  and  SZM  use  an  Asselin  filter  to  suppress  the  time  splitting  that  can  occur  with  leapfrog 
(Asselin  1972).  The  Asselin  filter  is  applied  at  the  end  of  each  timestep  to  the  fields  at  time  level 
(n)  by  averaging  in  a  bit  of  the  values  at  the  (n  +  1)  and  ( n  -  1)  time  levels,  i.e., 


$<*)  =  v(4>(*  +1)  +  <|><"  - 1>)  +  (1  -  2v)(j)(")  ,  (28) 

where  v  is  the  filter  coefficient.  If  (28)  is  rewritten  as 


<()(«)  =  <(>(«)  +  v(<j)(n  +1>  -  2<J)(">  +  4><n  - x)) , 


(29) 
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the  filter  looks  like  a  numerical  diffusion  term,  which  is  the  way  that  it  behaves.  A  typical  value 
of  v  used  in  POM  and  SZM  is  0.05,  and  this  is  the  value  used  in  the  tests  in  this  report. 

ECOM-si  uses  a  two-time-level  scheme  in  which  most  of  the  terms  in  the  baroclinic  equations 
are  evaluated  at  the  current  (n)  time  level.  This  is  sometimes  referred  to  as  a  forward  time-differencing 
scheme.  Vertical  diffusion,  however,  is  treated  fully  implicitly  as  in  POM  and  SZM.  The  temporal 
differencing  of  the  u  momentum  equation  for  ECOM-si  can  be  written  as 


- - — —  =  -  V  •  (v«)<»>  +  /V<">  -  +F(n)  +  - 

At  v  '  J  p n  dx  “  dz 


K, 


du(n+V 


AT 


dz 


(30) 


Note  that  ECOM-si  updates  the  velocities  before  updating  T  and  5  and  uses  the  new  velocities  to 
advect  T  and  S  so  that  in  the  T  and  5  equations,  the  advection  velocity  v  is  at  (n  +  1)  rather  than 
at  ( n )  as  in  the  momentum  equations. 

There  are  some  advantages  to  the  forward  scheme:  (a)  prognostic  variables  only  need  to  be 
saved  at  a  single  time  level,  which  saves  computer  storage,  (b)  there  is  no  time-splitting  problem 
that  one  must  deal  with  as  there  is  when  using  the  leapfrog  scheme,  and  (c)  timestep  restrictions 
are  calculated  for  a  timestep  At,  instead  of  with  2 At,  which  is  the  effective  timestep  for  the  leapfrog 
scheme.  The  big  disadvantage  of  the  forward  scheme  is  that  the  terms  are  not  centered  in  time  and 
can  suffer  significant  temporal  truncation  error.  The  magnitude  of  this  truncation  error  depends  on 
the  ratio  of  the  model  timestep  to  the  timescale  of  the  physical  process  being  simulated.  The  advection 
terms  for  the  forward  scheme  require  a  certain  amount  of  diffusion  to  suppress  numerical  dispersion 
and  can  become  very  dispersive  when  Ata2/(2A*)  >  1  (Clancy  1981),  where  A*  is  the  horizontal 
diffusion  coefficient.  This  is  discussed  further  in  the  tests  of  the  advection  schemes  in  Sec.  3.0. 


SCRUM  calculates  the  temporal  derivatives  for  the  baroclinic  equations  using  the  difference 
between  successive  time  levels,  i.e.,  the  value  at  ( n  +  1)  minus  the  value  at  (n),  as  does  ECOM-si. 
However,  SCRUM  uses  a  third-order  Adams  Bashforth  (AB3)  scheme  for  the  baroclinic  pressure 
gradient,  horizontal  advection,  and  Coriolis  terms.  This  scheme  provides  third-order  temporal  accuracy 
by  using  a  quadratic  extrapolation  of  the  values  at  the  three  previous  time  levels  to  the  value  at 
(n  +1),  so  that  the  calculated  value  is  centered  with  respect  to  the  time  derivative.  A  Crank-Nicolson 
scheme  is  used  for  vertical  advection  and  vertical  diffusion,  where  half  the  term  is  evaluated  at 
( n  +  1)  and  half  at  (n).  The  temporal  differencing  of  the  w-momentum  equation  for  SCRUM  can  be 
written  as 


u^*  1)-«(") 
At 


,  „  .  (a)  ,  (a)  1  dP  _(«)  d 

“  dz  1  P  n  dx  «  dz 


(«) 


K; 


du& 


dz 


(31) 


where  the  superscript  (a)  refers  to  the  AB3  scheme  where 

„<«)  =  !!„<»)  «  „(»-!)  +  5  „<»-2) 

12  12  12 

and  (c)  refers  to  the  Crank-Nicolson  scheme  where 


(32) 


(33) 
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ECOM-si’s  forward  treatment  of  the  Coriolis  term  is  numerically  unstable  in  that  inertial 
oscillations  grow  exponentially  with  an  e-folding  timescale  of  2/(A tf^).  If  the  timestep  used  in 
ECOM-si  is  small,  say  less  than  about  200  s,  which  is  generally  the  case  for  high-resolution  coastal 
simulations,  the  timescale  for  the  growth  of  the  inertial  oscillations  is  fairly  slow,  i.e.,  more  than 
10  d.  Such  a  slow  growth  of  inertial  motions  may  be  masked  by  other  variability  in  the  model  or 
cancelled  by  the  natural  damping  of  the  inertial  motions  due  to  friction  and  dispersion.  For  larger 
timesteps,  however,  the  growth  of  the  inertial  motions  can  be  a  significant  problem,  e.g.,  for 
A t  =  1200  s,  the  timescale  is  3.6  d  at  30°  N.  For  the  purpose  of  conducting  this  study,  which 
involved  performing  some  simulations  on  fairly  coarse  grids,  the  Coriolis  terms  in  ECOM-si  were 
converted  to  a  second-order  Adams-Bashforth  (AB2)  treatment.  For  the  u  equation,  this  is  given  by 


„(«  +  !)_„(«) 
A  t 


(34) 


The  AB2  scheme  provides  a  linear  extrapolation  of  the  Coriolis  terms  to  ( n  +i)  so  that  they  are 
centered  in  time  with  respect  to  the  time  derivative.  With  this  treatment  of  the  Coriolis  terms,  the 
amplification  of  inertial  oscillations  is  neutral,  i.e.,  there  is  no  significant  growth  or  damping. 
From  the  computational  point  of  view,  the  cost  of  the  AB2  treatment  of  the  Coriolis  terms  is  very 
small.  This  modification  to  ECOM-si  was  used  for  all  the  tests  conducted  in  this  report  and  no 
drawbacks  to  this  change  were  observed.  In  several  tests,  comparison  with  the  original  treatment 
of  the  Coriolis  terms  verified  the  improvement  in  stability  and  accuracy  with  the  AB2  treatment. 
Comparison  of  the  AB2  scheme  with  an  AB3  scheme  for  the  Coriolis  terms  in  some  of  the  tests 
showed  similar  results  for  these  two  schemes. 


2.7  Numerical  Treatment  of  the  Free  Surface  Mode 

A  limitation  of  most  explicit  numerical  schemes  is  that  signals  cannot  propagate  more  than  a 
single  grid  interval  in  a  single  timestep,  i.e.,  the  timestep  is  restricted  by 


(35) 


where  c  is  the  speed  of  the  signal  (for  leapfrog  schemes,  the  timestep  is  effectively  2A t  so  the 
restriction  is  At  <  Ax/(2 c)). 


Surface  gravity  waves  travel  very  quickly,  e.g.,  in  water  of  100-m  depth,  the  speed  of  surface 

gravity  waves  is  (g//)2  =  32  m/s.  For  a  grid  spacing  of  Ax  =  1  km,  explicit  treatment  of  these  waves 
in  an  ocean  model  requires  a  timestep  of  less  than  30  s.  If  the  surface  gravity  waves  are  excluded, 
the  timestep  is  limited  by  the  propagation  of  internal  waves  and  advection  velocities.  Since  these 
speeds  are  generally  less  than  about  2  m/s,  the  timestep  limitation  for  internal  wave  propagation  and 
advection  for  the  same  1-km  grid  would  be  about  500  s  (250  s  for  a  leapfrog  scheme). 

The  timestep  restriction  for  explicit  treatment  of  the  surface  waves  ranges  from  10  to  more 
than  50  times  smaller  than  the  timestep  restriction  for  the  other  processes  in  the  model,  depending 
on  the  water  depth,  the  maximum  internal  wave  speed,  and  the  maximum  current  velocities.  Hence, 
the  speed  of  surface  gravity  waves  will  severely  limit  the  timestep  of  an  ocean  model  unless  a 
special  numerical  treatment  is  applied  to  these  waves. 
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The  models  here  use  two  approaches  to  avoid  the  timestep  restriction  of  surface  gravity  waves. 
POM  and  SCRUM  use  the  so-called  “split-explicit”  scheme,  wherein  the  terms  governing  the 
propagation  of  the  surface  gravity  waves  are  treated  explicitly,  but  with  a  much  smaller  timestep 
than  the  rest  of  the  model  (Blumberg  and  Mellor  1987;  Song  and  Haidvogel  1994;  Killworth  et  al. 
1991;  Mellor  1996).  With  the  split-explicit  scheme,  the  3-D  momentum  and  continuity  equations 
are  vertically  averaged  to  obtain  the  barotropic,  2-D  momentum  and  continuity  equations  that 
govern  the  propagation  of  the  surface  gravity  waves.  The  vertically  averaged  equations  are  solved 
explicitly  with  a  sufficiently  small  timestep  to  avoid  violating  the  timestep  restriction  for  the 
surface  gravity  waves.  After  the  new  surface  elevation  and  depth-averaged  velocities  have  been 
solved  by  taking  a  number  of  small  timesteps,  the  3-D  equations  are  solved  with  a  much  larger 
timestep  governed  by  the  timestep  restriction  for  internal  waves  and/or  advection.  POM  and  SCRUM 
use  leapfrog  time  differencing  for  the  barotropic  mode.  POM  uses  an  Asselin  temporal  filter  to 
suppress  time  splitting  of  the  barotropic  equations  and  SCRUM  uses  a  trapezoidal  correction. 

ECOM-si  and  SZM  avoid  the  timestep  restriction  of  the  surface  gravity  waves  by  treating  the 
primary  terms  governing  the  propagation  of  these  waves  implicitly.  These  are  the  surface  pressure 
gradient  terms  in  the  momentum  equations  and  the  transport  gradient  terms  in  the  depth-averaged 
continuity  equation.  In  ECOM-si,  these  terms  are  treated  fully  implicitly,  i.e.,  fully  at  the  new  time 
level.  In  SZM,  the  temporal  weighting  of  these  terms  can  be  set  by  the  user.  For  the  tests  conducted 
in  this  report,  an  equal  weighting  of  these  terms  at  the  old  (n  -  1)  and  new  (n  +  1)  time  levels  was 
used. 

A  disadvantage  of  the  implicit  scheme  with  respect  to  the  split-explicit  scheme  for  the  surface 
gravity  waves  is  that  the  implicit  scheme  has  much  less  temporal  accuracy  since  it  is  using  a  much 
larger  timestep.  The  loss  of  accuracy  is  most  severe  at  smaller  wavelengths.  An  advantage  of  the 
implicit  scheme  is  that  it  tends  to  be  more  computationally  efficient,  especially  if  part  of  the  model 
domain  is  very  deep,  so  that  the  timestep  for  the  explicit  treatment  of  the  surface  waves  must  be 
very  small.  Another  advantage  of  the  implicit  scheme  over  the  split-explicit  scheme  is  that  it  is 
simpler  to  implement.  The  split-explicit  scheme  must  be  implemented  carefully  to  eliminate  incon¬ 
sistencies  between  the  depth-averaged  equations  and  the  rest  of  the  model.  Potential  problems  are 
discussed  by  Killworth  et  al.  (1991),  Dukowicz  and  Smith  (1994),  and  Mellor  (1996). 


2.8  Parameterization  of  Advection 

POM,  ECOM-si,  and  SZM  use  the  flux  form  of  the  advection  terms,  as  written  in  (1),  (2),  (6), 
and  (7).  SCRUM  uses  the  advective  form  in  which  the  advection  term  for,  e.g.,  temperature,  is 
written  as  v  •  VT  rather  than  V  •  (vT). 

An  advantage  of  implementing  the  advection  terms  in  flux  form  is  that  the  advected  field  can 
be  conserved.  However,  to  ensure  strict  conservation,  it  is  essential  that  the  numerical  form  of  the 
time  derivative  and  the  advection  terms  for  the  field  being  advected  be  consistent  with  the  numerical 
form  of  the  continuity  equation.  A  simple  way  to  check  for  this  in  the  models  is  to  set  the  advected 
field  to  a  constant.  Then  the  time  derivative  and  advection  terms  for  the  advected  field  must  cancel 
numerically  or  there  will  be  spurious  sources  and  sinks  of  the  advected  field  that  will  occur  at  grid 
cells  in  which  the  net  advective  volume  flux  into  a  grid  cell  and  the  change  in  volume  of  the 
grid  cell  do  not  match. 

POM,  ECOM-si,  and  SZM  go  to  some  trouble  to  ensure  that  the  advection  terms  for  T  and  S 
are  consistent  with  the  numerical  form  of  the  continuity  equation.  However,  this  is  not  quite  the 
case  for  advection  of  momentum.  To  have  fully  conservative  momentum  advection  would  require 
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an  iterative  procedure,  since  the  momentum  advection  terms  are  needed  to  calculate  the  new 
surface  elevation  and  the  new  surface  elevation  is  needed  to  calculate  the  momentum  advection 
terms  in  a  completely  conservative  manner.  (Another  smaller  source  of  error  in  calculating  momentum 
advection  is  that,  with  a  curvilinear  grid  and  grid-cell  thicknesses  that  change  with  time,  interpolation 
of  nondivergent  advective  volume  fluxes  to  the  staggered  momentum  grid  cells  on  a  C  grid  will 
not  result  in  exact  nondivergence  at  the  momentum  grid  cells.)  SZM  provides  for  such  an  iteration. 
However,  use  of  this  option,  which  adds  considerably  to  the  model’s  run  time,  does  not  generally 
have  a  significant  effect  on  the  model  results.  This  provides  some  validation  of  the  assumption 
made  in  these  models  that  strict  conservation  of  momentum  is  generally  less  critical  to  model 
accuracy  than  conservation  of  scalar  fields  because  of  the  less  dominant  role  of  advection  in  the 
momentum  equations. 

With  the  advective  form  of  the  model  equations  used  by  SCRUM,  strict  numerical  conservation 
of  the  advection  field  is  not  required  to  avoid  spurious  numerical  effects.  The  disadvantage  of  the 
advective  form  is  that  the  field  being  advected  is  not  strictly  conserved.  Note  that  the  strict 
conservation  of  scalar  fields  is  desirable,  but  tends  to  be  more  of  a  concern  for  long  time  integrations 
than  the  typical  short  simulations  usually  conducted  in  coastal  modeling,  and  also  more  of  a 
concern  for  large  ocean  regions  than  for  coastal  areas  where  sources,  sinks,  and  flows  through  open 
boundaries  play  a  larger  role.  Of  course,  errors  due  to  nonconservation  of  advected  fields  must 
remain  sufficiently  small  so  as  not  to  significantly  affect  results. 


2.9  Parameterization  of  Horizontal  Mixing 

POM,  ECOM-si,  and  SZM  use  Laplacian  horizontal  mixing.  POM  and  ECOM-si  use  the 
Smagorinsky  (1963)  parameterization  in  which  the  horizontal  friction  terms  are  written  as 
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For  the  Smagorinsky  mixing  scheme,  the  horizontal  mixing  coefficient  AM  is  calculated  as  a 
function  of  the  local  horizontal  grid  resolution  and  velocity  shear,  i.e., 


AM=CsAxAy 
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(38) 


The  magnitude  of  the  Smagorinsky  eddy  coefficients  is  scaled  by  the  constant  Cs  (called  HORCON 
in  the  POM  and  ECOM-si  model  programs).  Values  used  for  Cs  range  from  about  0.01  to  0.5. 
Large  values  of  Cs  tend  to  dissipate  small  features,  whereas  values  that  are  too  small  can  lead  to 
numerical  noise  and/or  instability.  Typically,  a  value  near  0.1  is  used,  and  0.1  is  the  value  used 
for  the  tests  conducted  in  this  report. 
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The  horizontal  mixing  coefficient  Ay  in  SZM  is  set  equal  to  the  maximum  of  a  constant  background 
value  Aa  and  a  value  needed  to  keep  the  grid-cell  Reynolds  Number  Re  below  a  maximum 
specified  value,  i.e.,  4 


(u  Ax 


(41) 


The  value  of  Re  is  typically  set  to  a  value  of  30-100  in  the  horizontal.  Values  in  this  range 
generally  provide  a  good  compromise  between  allowing  natural  dynamic  instability  processes  to 
develop  and  keeping  numerical  noise  to  acceptable  levels.  In  most  of  the  tests  conducted  here,  Re 
is  set  to  a  value  of  100  and  Aa  =  aQ  Ax,  where  Aa  =  0.1  cm/s. 

The  form  used  for  horizontal  diffusion  of  T  and  5  in  POM,  ECOM-si,  and  SZM  is  the  Laplacian 
form  expressed  in  (39)  and  (40).  The  models  provide  for  allowing  Afj  to  be  some  fraction  of  AM, 
but  are  frequently  run  with  Afj  =  Aj^.  Ajj  and  AM  were  set  equal  in  the  tests  conducted  here. 

In  regions  where  there  are  large  changes  in  the  bottom  depth,  the  horizontal  diffusion  along 
sloping  sigma  layers  can  cause  severe  cross-isopycnal  diffusion  (Paul  1994).  An  expedient  that  is 
sometimes  used  with  sigma  coordinate  models,  which  is  employed  in  the  models  tested  here,  is  to 
subtract  a  smooth,  background  field  from  the  T  ox  S  fields  when  calculating  horizontal  diffusion. 
By  calculating  the  horizontal  diffusion  based  on  the  anomaly  from  a  smooth  background  field,  most 
of  the  component  of  vertical  diffusion  that  occurs  when  diffusion  is  calculated  along  sloping  sigma 
layers  is  eliminated.  In  the  tests  conducted  here,  the  background  T  and  5  fields  are  calculated  as 
the  horizontally  averaged  T  and  S  profiles.  This  is  a  commonly  used  procedure.  An  alternative 
strategy  is  to  use  smooth  but  horizontally  varying  background  fields  to  accommodate  changes  in  the 
structure  of  the  T  and  S  fields  in  different  parts  of  the  model  domain.  The  background  fields  can 
also  be  periodically  updated  to  accommodate  changes  in  T  and  S  that  occur  in  time.  We  note  that 
the  use  of  these  procedures  can  significantly  reduce  the  problem  of  severe  cross-isopycnal  diffusion; 
however,  they  can  introduce  other  problems,  one  example  of  which  will  be  discussed  in  the 
upwelling/downwelling  test  in  Sec.  8.0. 

SCRUM  offers  both  Laplacian  horizontal  friction  (with  constant  values  of  AM  and  AH)  and 
high-order,  biharmonic  friction  in  the  horizontal.  Biharmonic  friction  is  more  effective  in  removing 
small-scale  noise  than  Laplacian  friction,  but  is  not  strictly  conservative. 


2.10  Parameterization  of  Vertical  Mixing 

POM  and  ECOM-si  use  the  MYL2.5  turbulence  model  to  parameterize  vertical  mixing  (Mellor 
and  Yamada  1982;  Blumberg  and  Mellor  1987).  This  turbulence  model  requires  the  prognostic 
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calculation  of  two  turbulence  quantities,  q2,  which  is  twice  the  turbulent  kinetic  energy  (TKE),  and 
q2t,  where  (  is  the  turbulent  length  scale. 

SZM  uses  the  MYL2  turbulence  model  (Mellor  and  Yamada  1974;  Mellor  and  Durbin  1975). 
The  MYL2  turbulence  model  assumes  a  local  balance  between  shear  production,  buoyancy  production, 
and  dissipation  of  TKE,  and  does  not  require  prognostic  prediction  of  any  turbulence  quantities. 
This  saves  computer  time  and  memory  relative  to  using  the  MYL2.5  model;  however,  the  simplified 
TKE  equation  does  not  account  for  storage  and  transport  of  TKE.  The  turbulence  length  scale  t  for 
SZM  is  calculated  diagnostically  based  on  the  distance  from  the  top  ( zt )  and  bottom  (. zb )  of  a 
turbulent  boundary  layer,  i.e., 


1+zb1)  ’  (42) 

where  k  =  0.4  is  von  Karman’s  constant.  The  simplification  of  the  MYL2  scheme  relative  to  the 
MYL2.5  scheme  may  have  some  disadvantages  in  terms  of  accuracy  and  robustness.  Tests  will  be 
conducted  to  determine  how  the  implementations  of  the  MYL2  and  MYL2.5  schemes  in  the  models 
compare  for  some  basic  mixing  situations. 


2.11  Position  of  Variables  on  the  Model  Grid 

The  horizontal  positioning  of  the  main  variables  with  respect  to  the  grid  cells  is  the  same  for 
all  the  models.  They  all  use  an  Arakawa  C  grid  in  which  the  horizontal  velocities  are  defined  at 
the  edges  of  the  grid  cells,  and  the  other  variables  ( T ,  S,  p,  p,  w,  £)  are  defined  at  the  center  of  the 
grid  cells  (Fig.  1). 

Looking  at  the  vertical  location  of  the  variables  for  POM,  ECOM-si,  and  SZM,  all  the  main 
variables  are  defined  at  the  center  of  a  layer  except  for  w,  which  is  defined  at  the  interface  between 
the  layers  (Fig.  2a).  For  SCRUM,  which  uses  a  FE  scheme  in  the  vertical,  all  the  main  variables 
are  defined  at  the  layer  interfaces  (Fig.  2b).  For  POM  and  ECOM-si,  the  turbulence  variables  q2 
and  q2Q  are  defined  at  the  w  locations. 
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Fig.  1  —  Horizontal  arrangement  of  primary  model  variables  on 
the  numerical  grid  for  all  the  models 
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(b) 


Fig.  2  —  Vertical  arrangement  of  primary  model  variables  on  the  numerical  grid  for  (a)  POM, 

ECOM-si,  and  SZM  and  (b)  SCRUM 


The  horizontal  eddy  coefficients  for  POM  and  ECOM-si  are  defined  at  the  grid-cell  centers. 
For  SZM,  separate  horizontal  eddy  coefficients  are  calculated  for  diffusion  in  the  x  and  y  directions, 
and  these  are  defined  at  the  u  and  v  locations,  respectively.  The  vertical  eddy  coefficients  for  POM, 
ECOM-si,  and  SZM  are  defined  at  the  w  locations;  however,  for  SCRUM  they  are  defined  at  the 
center  of  the  grid  cells. 


2.12  Setting  Up  the  Models  for  Testing 

POM,  as  obtained  from  Princeton  University,  was  set  up  for  a  closed  domain,  and  the  do-loop 
limits  used  in  the  model  were  unsuitable  for  use  with  open  boundaries.  Hence,  before  beginning 
testing  of  POM,  all  the  do-loops  in  the  model  were  checked,  and  if  necessary,  modified  to  allow 
the  use  of  open  boundaries  on  all  sides  (note  that  the  do-loop  limits  needed  for  open  boundaries 
are  not  incompatible  with  the  use  of  either  closed  or  periodic  boundaries). 

This  problem  of  the  do-loop  limits  has  been  noted  by  many  users  of  POM.  In  1996,  the  version 
of  POM  available  from  Princeton  was  replaced  with  a  new  version  in  which  the  do-loop  limits  are 
reported  to  have  been  modified  to  more  easily  allow  the  use  of  open  boundaries  (this  new  version 
of  POM  has  not  yet  been  compared  with  our  own  modifications). 

ECOM-si  came  with  a  fairly  well-developed  user  interface  that  allows  setting  up  a  model  run 
for  a  particular  simulation  by  placing  all  the  input  parameters  and  data  required  for  the  run  in  two 
input  data  files.  An  advantage  of  this  setup  is  that  few  (if  any)  changes  need  to  be  made  to  the 
model  code  to  run  a  particular  simulation.  However,  this  setup  proved  to  be  cumbersome  for 
the  tests  conducted  here,  and  was  replaced  with  a  set  of  subroutines  to  define  the  model  setup, 
initialization,  and  forcing. 

The  ability  to  use  periodic  lateral  boundary  conditions  in  both  x  and  y  was  added  to  both  POM 
and  ECOM-si  (this  feature  was  already  available  for  SZM  and  SCRUM).  Periodic  boundaries  are 
not  needed  for  simulations  of  actual  coastal  regions,  but  are  very  useful  for  testing  purposes.  Many 
of  the  tests  conducted  for  this  model  comparison  study  used  periodic  boundaries  in  one  or  both 
directions. 

As  discussed  in  the  introduction,  there  were  a  number  of  problems  found  with  the  Version  2.1 
SCRUM  code  that  was  obtained  from  Rutgers.  The  Crank-Nicolson  treatment  of  vertical  mixing 
tended  to  produce  noisy,  inaccurate  vertical  mixing,  and  some  terms  were  missing  from  the  forcing 
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for  the  free  surface  mode.  These  problems  were  generally  due  to  the  fact  that  SCRUM  was  still 
undergoing  changes.  As  previously  stated,  because  of  these  problems  and  extensive  changes  that 
were  made  to  SCRUM  in  the  more  recent  version  of  the  model  that  has  been  made  available 
(Version  3.0),  SCRUM  was  not  included  in  most  of  the  model  tests. 


3.0  TEST  OF  HORIZONTAL  ADVECTION 
3.1  Rotating  Cone  Advection  Test 

The  performance  of  the  horizontal  advection  schemes  used  in  the  models  was  investigated  by 
performing  the  classic  “rotating  cone”  advection  experiment  (Crowley  1968;  Molenkamp  1968; 
Orszag  1971;  Anderson  and  Gattahi  1974;  Purnell  1976;  Shannon  1979;  Long  and  Pepper  1981; 
McRae  et  al.  1982;  Smolarkiweicz  and  Grabowski  1990;  Cantekin  and  Westerink  1990).  In  this 
experiment,  a  cone-  or  Gaussian-shaped  feature  is  advected  in  a  circular  path  by  an  advection  field 
representing  solid  body  rotation  about  the  center  of  the  domain.  A  perfect  advection  scheme  would 
advect  the  feature  in  a  circular  path  at  the  proper  speed  without  changing  its  shape.  The  performance 
of  the  advection  scheme  is  judged  by  the  distortion  of  the  shape  of  the  feature  after  a  certain  period 
of  advection  and  by  errors  in  its  position. 

Figure  3a  shows  the  model  domain  and  the  initial  condition  for  this  experiment.  The  model 
domain  consists  of  a  square  region  40  km  on  a  side.  The  grid  spacing  within  the  domain  is  a 
uniform  1  km.  The  initial  feature,  which  is  taken  to  be  in  the  shape  of  a  cone,  is  located  10  km 
south  of  the  center  of  the  model  domain  and  has  a  radius  at  its  base  of  4  km  (4  gridpoints)  and 
a  maximum  amplitude  of  50  units.  The  field  being  advected  could  be  considered  to  be  temperature, 
salinity,  or  some  other  constituent  of  the  model.  The  advection  field  consists  of  solid  body  rotation 
about  the  center  of  the  model  domain  with  a  period  of  rotation  of  1.6  d.  With  this  rotation  period, 
the  advection  velocity  at  the  location  of  the  center  of  the  cone  is  45.4  cm/s. 

Note  that  the  advection  tests  were  conducted  with  the  numerical  schemes  used  for  advection 
in  the  models,  but  not  with  the  models  themselves.  This  was  the  only  test  in  this  report  that  was 
conducted  in  this  manner.  The  advection  test  was  conducted  this  way  because  it  was  simpler  to  do 
so  than  to  set  up  this  test  with  the  models. 

These  advection  experiments  included  a  certain  amount  of  explicit  diffusion.  This  was 
necessary  since  some  diffusion  is  needed  to  filter  the  numerical  noise  generated  by  the  advection 
schemes.  For  the  advection  experiments,  the  diffusion  terms  were  parameterized  using  standard, 
second-order,  flux-conserving,  Laplacian  diffusion  with  a  constant  diffusion  coefficient. 

As  discussed  in  the  previous  section,  POM,  ECOM-si,  and  SZM  use  the  same  second-order, 
centered,  spatial  finite  difference  scheme  for  the  horizontal  advection  terms.  With  this  scheme,  the 
field  being  advected  is  interpolated  from  the  grid-cell  centers  to  the  grid-cell  boundaries  using 
two-point  averages,  and  the  flux  through  the  boundaries  of  the  cell  is  computed  by  multiplying  the 
value  of  the  field  at  the  grid-cell  boundary  by  the  normal  velocity  at  the  boundary.  The  change  in 
the  field  due  to  advection  is  then  calculated  as  the  sum  of  the  fluxes  into  and  out  of  each  of  the 
grid  cells.  This  scheme  is  conservative,  since  the  flux  leaving  a  cell  is  equal  to  the  flux  entering 
an  adjoining  cell  though  the  common  face. 

SCRUM  uses  the  advective  form  of  the  advection  term  in  which  the  gradient  of  the  field 
being  advected  is  calculated  at  the  grid-cell  center  using  a  two-point  difference  and  is  multiplied 
by  the  advection  velocity  averaged  to  the  grid-cell  center  using  a  two-point  average  (since  the  advection 


Fig.  3  —  Results  of  rotating  cone  advection  test  after  one  revolution  (1.6  d)  for  (b)  leapfrog  scheme  used  by 
POM  and  SZM,  (c)  AB3  scheme  used  by  SCRUM,  and  (d)  forward  temporal  scheme  used  by  ECOM-si.  The 
initial  position  and  shape  of  the  cone  are  shown  in  (a).  The  initial  height  of  the  cone  is  50  units.  The  timestep 
was  360  s  for  the  leapfrog  and  AB3  schemes  and  180  s  for  the  forward  scheme.  The  horizontal  viscosity  was 
a  constant  10  m2/s.  The  gridpoint  locations  are  denoted  by  the  tick  marks,  and  Ajc  =  1  km. 


field  is  constant  in  time  and  varies  gradually  in  space,  the  flux  and  advective  forms  of  the 
advection  term  gave  very  similar  results  in  this  test). 

The  differences  in  the  behavior  of  the  advection  schemes  used  by  these  models  in  this  test  are 
due  mainly  to  the  differences  in  the  temporal  differencing  schemes.  As  discussed  in  Sec.  2.6,  POM 
and  SZM  use  a  leapfrog  temporal  scheme  in  which  the  advection  terms  are  calculated  at  the  central 
time  level,  i.e,  for  the  advection  of  the  variable  J,  and  ignoring  other  terms,  we  have 


Tin  +  1)  =  T(n  - 1)  +  2A teW  , 


(43) 
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where  denotes  the  advection  term.  The  temporal  accuracy  of  leapfrog  advection  is  second 
order,  i.e.,  the  temporal  truncation  errors  are  proportional  to  (At)2. 

ECOM-si  uses  the  two-time-level,  forward  temporal  scheme  in  which  the  advection  terms  are 
evaluated  at  ( n ),  i.e., 


r(«  +  !)  =  j(")  +  AtCW  .  (44) 

The  temporal  accuracy  of  the  forward  scheme  is  only  first  order  because  the  advection  terms  are 
not  centered  in  time. 

SCRUM  updates  the  baroclinic  fields  from  the  values  at  (n)  using  the  AB3  scheme 

T(n  + 1)  _  T(n)  +  At((23 CW  -  16C(n~V  +  5C(n-2))/12)  ,  (45) 

for  which  the  truncation  error  is  proportional  to  (At)3. 

The  standard  timestep  for  the  advection  experiments  was  taken  to  be  At  =  360  s  and  the 
horizontal  diffusivity  AH  was  set  to  a  constant  10  m2/s.  A  requirement  for  numerical  stability  of 
an  explicit  advection  scheme  is  that 


_  At 
Cu~UAx<1, 


(46) 


where  Cu  is  the  Courant  Number.  A  similar  stability  requirement  for  explicit  diffusion  is  given  by 
the  diffusion  number  d  (Roache  1976). 


d  =  2AH-^<l.  (47) 

Axa 

For  the  leapfrog  scheme,  the  value  of  At  in  (46)  and  (47)  should  be  replaced  by  2At,  since  the 
timestep  is  effectively  2A t. 

Also,  as  mentioned  in  Sec.  2.0,  it  is  frequently  desirable  for  the  grid-cell  Reynolds  Number 
(Re  =  uAx/Ah),  which  is  the  ratio  of  the  relative  strengths  of  the  advection  and  diffusion  terms, 
to  be  in  the  range  of  30-100  so  that  the  natural  physical  instabilities  of  the  flow  on  the  smallest 
spatial  scales  resolvable  by  the  model  are  not  damped  by  the  diffusion  terms. 

For  the  temporal  differencing  schemes  used  by  ECOM-si  and  SCRUM,  and  for  our  standard 
parameter  values  (At  =  360  s,  Ax  =  1  km,  and  Afj  =  10  m2/s),  Cu  =  0.16,  d  =  0.0072,  and  RSg  =  45 
at  the  location  of  the  center  of  the  cone  where  u  =  45.4  cm/s.  In  the  corners  of  the  model  grid 
where  u  -  128.5  cm/s,  Cu  =  0.46  and  Reg  =  127.  For  the  leapfrog  scheme  used  by  POM  and  SZM, 
the  timestep  is  effectively  2  x  360  s  in  (46)  and  (47)  so  that  values  of  Cu  are  twice  as  large,  i.e.,  0.33  at 
the  location  of  the  center  of  the  cone  and  0.93  in  the  comers.  Hence,  the  schemes  satisfy  the  stability 
constraints  described  by  (46)  and  (47)  over  the  entire  grid  for  the  standard  parameter  values. 

For  the  forward  advection  scheme  used  by  ECOM-si,  there  is  another  constraint  (Clancy 
1981),  which  is 
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A  tu‘ 


<2 


lH 


(48) 


Failure  to  comply  with  this  restriction  results  in  severe  numerical  dispersion  and  possible  instability. 
Note  that  for  the  standard  parameters  values,  Atu2/AH=7.3,  i.e.,  (48)  is  not  satisfied.  Equation  (48) 
could  be  satisfied  by  reducing  the  ratio  A t/AH  by  a  factor  of  about  4. 

Figure  3  shows  the  initial  shape  of  the  cone  and  the  results  for  the  advection  of  the  cone  after 
one  complete  revolution  (1.6  d)  for  each  of  the  three  advection  schemes:  the  leapfrog  advection 
scheme  used  by  POM  and  SZM,  the  AB3  scheme  used  by  SCRUM,  and  the  forward  scheme  used 
by  ECOM-si.  Results  from  all  the  advection  experiments  are  also  summarized  in  Table  2.  The 
forward  advection  scheme  was  unstable  with  At  =  360  s  and  AH  =  10  m2/s,  so  the  timestep  was 
reduced  by  half,  i.e.,  the  results  for  the  forward  scheme  in  Fig.  3  are  for  At  =  180  s. 

The  results  for  the  leapfrog  and  AB3  schemes  in  Fig.  3  are  about  the  same.  The  peak  of  the 
cone  is  reduced  from  50  to  19  units  due  to  a  combination  of  dispersion  caused  by  the  numerical 
scheme  (i.e.,  failure  of  the  scheme  to  advect  the  field  at  the  true  advection  velocity),  and  actual 
diffusion  by  the  diffusion  term  (the  diffusion  term  alone  would  have  reduced  the  peak  from  50  to 
25  units  in  1.6  d).  Numerical  dispersion  creates  a  wake  behind  the  cone  that  has  a  maximum  value 
of  -5  units  in  the  first  wave  behind  the  cone.  Dispersion  is  also  responsible  for  the  fact  that  the 
peak  of  the  cone  does  not  quite  get  back  to  its  initial  position,  but  falls  short  by  about  3  km. 

Figure  4  shows  results  for  the  leapfrog  and  AB3  schemes  for  the  same  experiment,  but  with 
Ah  reduced  from  10  to  5  m2/s.  This  decrease  in  AH  gives  RSg  =  90  at  the  location  of  the  cone.  The 
leapfrog  and  AB3  schemes  again  give  similar  results.  With  the  reduced  diffusivity,  the  maximum 


Table  2  —  Results  from  Cone  Experiment  with  Advection  Schemes 
Used  by  Ocean  Models 


SCHEME 

At 

(s) 

DIFFUSIVITY 

(m2/s) 

MAX 

AMPLITUDE 

WITHIN 

CONE 

MAX 

AMPLITUDE 

WITHIN 

WAKE 

PHASE 

ERROR 

(km) 

Leapfrog 

360 

10 

19 

-5 

3 

Leapfrog 

360 

5 

22 

-8 

3 

Leapfrog 

90 

5 

22 

-8 

3 

AB3 

360 

19 

-5 

3 

AB3 

360 

5 

22 

-7 

3 

AB3 

90 

5 

22 

-8 

3 

Forward 

180 

10 

24 

-13 

3 

Forward 

90 

10 

21 

-8 

3 

Forward 

45 

10 

20 

-6 

3 

Forward 

30 

10 

19 

-6 

3 

Forward 

20 

10 

19 

-6 

3 

Forward 

180 

20 

18 

-7 

3 
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Fig.  4  —  Results  of  rotating  cone  advection  test  after  one  revolution  (1.6  d)  for  (a)  leapfrog  and  (b)  AB3 
schemes  with  a  constant  horizontal  viscosity  of  5  m2/s 


amplitude  within  the  cone  is  increased  from  19  to  22  units,  but  the  amplitude  of  the  wake,  caused 
by  numerical  dispersion,  is  also  increased.  Reducing  the  timestep  with  these  two  schemes  did  not 
significantly  affect  the  results  (Table  2),  i.e.,  the  error  is  primarily  due  to  spatial,  rather  than 
temporal,  truncation  error. 

The  forward-in-time  advection  scheme  is  well  known  to  be  highly  dispersive  (Roache  1976) 
and  the  results  with  this  scheme  in  Fig.  3  show  much  greater  dispersion  than  with  the  leapfrog  and 
AB3  schemes,  in  spite  of  the  fact  that  a  smaller  timestep  is  used.  The  maximum  value  in  the  wake 
behind  the  cone  after  one  revolution  is  -13  units,  compared  with  -5  units  for  the  leapfrog  and  AB3 
schemes. 

Figure  5a-c  shows  results  for  the  forward  scheme  with  timesteps  of  90,  45,  and  20  s  (with 
Aff  =  10  m2/s).  The  values  of  Cu  at  the  center  of  the  cone  for  these  experiments  are  0.04,  0.02,  and 
0.01,  i.e.,  well  below  the  value  of  1.0  required  for  stability  with  the  leapfrog  and  AB3  schemes. 
Because  of  the  large  temporal  truncation  error  of  the  forward  scheme,  the  dispersion  is  significantly 
reduced  as  the  timestep  is  reduced.  The  results  with  At  =  45  s  are  similar  to  the  results  with  the 
other  two  schemes.  Further  reduction  of  the  timestep  below  20  s  does  not  improve  the  results, 
indicating  that  the  remaining  error  is  primarily  due  to  spatial  truncation  error. 

Increasing  the  diffusivity  can  also  reduce  the  dispersion  of  the  forward  advection  scheme. 
Figure  5d  shows  a  result  with  At  =  180  s  as  was  used  for  the  result  in  Fig.  3d,  but  with  AH 
increased  to  from  10  to  20  m2/s.  Note  that  with  this  larger  value  of  Ajj,  the  forward  scheme 
satisfies  the  criteria  expressed  in  (48).  Figure  5d  illustrates  that  increasing  the  diffusivity  can  reduce 
the  numerical  dispersion,  though  at  the  cost  of  increasing  the  diffusiveness  of  the  numerical  scheme. 

In  summary,  the  advection  schemes  used  by  POM,  SZM,  and  SCRUM  give  similar  results  in 
these  tests.  The  numerical  dispersion  of  these  schemes  is  certainly  significant.  However,  based  on 
the  wide  use  of  these  advection  schemes  in  ocean  modeling,  this  degree  of  dispersion  has  been 
deemed  to  be  tolerable  (in  recent  years,  ocean  modelers  have  been  increasingly  investigating  the 
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Fig.  5  —  Results  of  rotating  cone  advection  test  after  one  revolution  (1.6  d)  with  forward  temporal 
scheme  for  different  timestep  and  viscosity  values 
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use  of  alternative  advection  schemes  (Hecht  et  al.  1995)).  The  error  of  the  advection  schemes  used 
in  POM,  SZM,  and  SCRUM  is  mostly  due  to  spatial  truncation  error  and  the  results  are  not  very 
sensitive  to  the  timestep  as  long  as  the  timestep  is  small  enough  to  keep  Cu  below  the  value 
required  for  numerical  stability. 

The  forward  advection  scheme  used  by  ECOM-si  is  well  known  to  be  highly  dispersive. 
Unlike  the  leapfrog  and  AB3  schemes,  the  forward  scheme  is  very  sensitive  to  the  timestep  because 
of  the  significant  temporal  truncation  error.  As  Cu  is  increased  above  a  value  of  about  0.04,  the 
dispersion  of  the  forward  scheme  increases  markedly.  This  dispersion  can  be  reduced  and  numerical 
stability  maintained  by  increasing  the  diffusivity  or  viscosity  so  that  A tu2/AH  remains  below  a 
value  of  about  two.  However,  as  the  value  of  Cu  increases  above  about  0.04,  the  values  of 
required  to  satisfy  this  criteria  will  result  in  increasingly  diffusive  flows. 
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Fig.  6  —  Results  of  advective  overshoot  experiment  with 
second-order,  centered  advection  scheme.  The  dashed  line 
shows  the  initial  temperature  and  the  solid  line  shows  the 
predicted  temperature  at  the  time  of  maximum  overshoot.  In 
(a),  the  initial  temperature  jump  occurs  between  two  gridpoints 
and  in  (b)  it  is  distributed  (linearly)  over  four  gridpoints.  The 
gridpoint  locations  are  indicated  by  tick  marks  along 
the  distance  axis  and  also  by  dots  at  the  predicted  temperature 
values. 


3.2  Advective  Overshoot  with  Second-Order,  Centered  Advection  Scheme 

The  second-order,  centered  advection  scheme  used  in  all  the  models  being  compared  here  is 
subject  to  an  error  sometimes  referred  to  as  “advective  overshoot,”  when  spatial  changes  in  the 
field  being  advected  are  not  adequately  resolved.  Figure  6  shows  the  results  of  the  uniform  advection 
of  a  field  in  which  there  is  a  sharp  change  (i.e.,  a  sharp  front)  occurring  between  gridpoints.  The 
field  being  advected  could  be  any  model  variable.  In  this  case,  the  advected  variable  is  taken  to 
be  temperature. 

In  the  case  shown  in  Fig.  6,  the  grid  spacing  is  1  km,  the  timestep  is  200  s,  and  the  advection 
velocity  is  a  constant  10  cm/s  in  a  direction  from  right  to  left.  The  value  of  the  horizontal  diffusivity 
has  been  taken  be  to  a  constant  1  m2/s  to  give  a  grid-cell  Reynolds  number  of  100.  The  dashed 
line  in  the  figure  shows  the  initial  value  of  the  advected  field  and  the  solid  line  shows  the  field 
after  70  timesteps.  In  Fig.  6a,  it  can  be  seen  that  the  value  at  the  gridpoint  on  the  upstream  side 
of  the  front  has  decreased  by  about  7°C  after  70  timesteps. 

The  reason  for  the  temperature  drop  on  the  upstream  side  of  the  front  is  illustrated  in  Fig.  7. 
With  the  second-order,  centered  advection  scheme,  the  value  of  the  field  advected  through  the  face 
of  a  grid  cell  is  the  mean  value  of  the  field  in  the  adjoining  grid  cells.  Hence,  the  value  of  the 
temperature  being  advected  out  of  the  cold  grid  cell  on  the  upstream  side  of  the  front  in  Fig.  7  is 
higher  than  the  value  being  advected  in  on  the  other  side,  and  the  value  within  the  grid-cell  drops. 
For  a  jump  in  the  advected  field  that  occurs  completely  between  two  grid  cells,  the  maximum  drop 
at  the  upstream  grid  cell  is  about  1/3  the  value  of  the  initial  jump.  The  size  of  this  drop  is  not 
significantly  affected  by  the  timestep  or  the  magnitude  of  the  velocity.  Nor  is  the  size  of  the  drop 
much  affected  by  the  horizontal  diffusivity,  as  long  as  the  diffusivity  is  sufficiently  small  that  the 
grid-cell  Reynolds  number  is  fairly  high  (which  is  generally  the  case  in  these  models  to  avoid 
excessive  diffusion),  e.g.,  if  the  horizontal  diffusivity  is  increased  by  a  factor  of  10  to  reduce  the 
grid-cell  Reynolds  number  to  10  (a  fairly  small  value  for  an  ocean  model),  the  drop  in  the  tem¬ 
perature  at  the  gridpoint  upstream  of  the  front  is  still  about  5°C.  The  results  of  several  experiments 
in  which  the  values  of  the  parameters  in  this  experiment  were  varied  are  summarized  in  Table  3. 
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TEMPERATURE  AT  THIS  CELL  DROPS 


Fig.  7  —  Schematic  showing  advection  between  grid  cells  at  the  location  of  a 
sharp  front.  The  temperature  flux  at  the  interface  between  the  grid  cells  is  calculated 
based  on  the  mean  temperature  of  the  adjoining  cells.  The  temperature  at  the  cold 
grid  cell  on  the  upwind  side  of  the  front  drops  because  the  advective  flux  leaving 
the  cell  is  much  warmer  than  the  advective  flux  entering  the  cell. 


Table  3  —  Results  of  Experiments  Conducted  to  Investigate  Advective  Overshoot  at 
a  Front  with  a  Second-Order,  Centered  Advection  Scheme.  The  Initial  (20°C)  Change 
in  the  Temperature  Field  Occurs  Completely  Between  Two  Grid  Cells.  The  Temperature 
Drop  Reported  in  the  Table  is  the  Maximum  Temperature  Drop  that  Occurs  at  the 
Cold  Grid  Cell  on  the  Upstream  Side  of  the  Front. 


GRID 

SPACING 

(km) 

TIMESTEP 

(s) 

VELOCITY 

(cm/s) 

DIFFUSIVITY 

(m2/s) 

CELL  RE 

MAX 

TEMPERATURE 
DROP  (°C) 

1 

100 

10.0 

100 

7.2 

1 

10 

1.0 

100 

7.2 

1 

1 

0.1 

100 

7.2 

1 

10 

10.0 

10 

4.9 

1 

Kfl 

10 

1.0 

100 

7.2 

The  magnitude  of  the  advective  overshoot  can  be  reduced  by  increasing  the  resolution  of  the 
front.  Figure  6b  shows  the  result  of  resolving  the  initial  front  by  four  gridpoints  instead  of  two. 
The  magnitude  of  the  overshoot  is  reduced  from  7.2  to  1.7°C. 

Sharp  temperature  or  salinity  fronts  between  adjacent  gridpoints  can  occur  in  model  simulations 
with  sigma  coordinates  when  there  is  a  sudden  change  in  depth  such  that  horizontally  adjacent 
gridpoints  in  the  bottom  layers  lie  on  different  sides  of  a  strong  thermocline  or  halocline.  Such  a 
situation  is  illustrated  in  Fig.  8.  In  this  case,  in  the  bottom  layer  of  the  model,  a  grid  cell  in  shallow 
water  on  the  shelf  is  adjacent  to  a  grid  cell  in  much  deeper  water  off  the  shelf  at  which  the  temperature 
is  20°C  colder.  If  there  is  advection  from  the  cold  grid  cell  to  the  warm  cell,  the  temperature  within 
the  cold  cell  will  drop  as  illustrated  in  Fig.  6a.  Such  a  situation  occurred  in  simulations  of  the  Gulf 
of  Mexico  circulation  with  POM  and  SZM  (Martin  1998).  The  resulting  advective  overshoot  caused 
the  bottom-layer  temperature  in  the  area  of  the  steep  continental  slope  south  of  Cuba  to  become 
negative  and  the  model  simulations  became  unstable.  The  solution  in  these  cases  is  to  increase  the 
horizontal  resolution  of  steep  bottom  slopes,  either  by  increasing  the  horizontal  grid  resolution  or 
by  reducing  the  slope  itself,  so  that  gradients  in  the  model  fields  are  better  resolved.  Note  that  this 
particular  problem  does  not  occur  with  z-level  vertical  coordinates,  since  the  gradient  between  the 
near  surface  and  deep  water  will  be  resolved  by  several  vertical  levels. 
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Fig.  8  —  Schematic  illustrating  the  large  change  in  temperature 
that  can  occur  between  adjoining  grid  cells  in  the  bottom 
layers  of  a  sigma  coordinate  model  if  there  is  a  large  change 
in  depth  between  the  cells  such  that  the  cells  lie  on  different 
sides  of  a  strong  thermocline 


4.0  TESTS  OF  VERTICAL  MIXING 

Tests  of  vertical  mixing  were  conducted  to  look  at  mixing  within  the  surface  and  bottom 
mixed  layers.  In  addition,  an  experiment  involving  the  formation  of  a  tidal  mixing  front  was 
performed  to  investigate  a  situation  where  the  surface  and  bottom  mixed  layers  combine  to  mix  the 
entire  water  column  from  top  to  bottom. 

The  POM,  ECOM-si,  and  SZM  models  use  a  fully  implicit  scheme  for  vertical  mixing.  With 
this  scheme,  the  vertical  diffusion  terms  are  evaluated  using  the  newly  calculated  values  of  the  field 
being  mixed.  For  example,  for  the  two-time-level  scheme  used  by  ECOM-si,  fully  implicit  mixing 
of  temperature  would  be  written  as 
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Since  the  new  ( n  +  1)  values  of  the  field  are  not  known  at  the  beginning  of  a  timestep,  the  implicit 
treatment  of  vertical  diffusion  couples  the  numerical  equations  in  the  vertical,  which  requires  the 
solution  of  a  tridiagonal  set  of  equations  at  each  horizontal  point. 

Although  the  fully  implicit  diffusion  scheme  suffers  a  loss  of  accuracy  from  temporal  truncation 
error  because  it  is  not  centered  in  time,  it  has  the  very  desirable  property  for  the  treatment  of 
vertical  diffusion  that  it  is  always  numerically  stable  and  physically  well-behaved.  As  a  practical 
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matter,  these  latter  properties  have  generally  been  considered  to  be  more  important  for  vertical 
diffusive  processes  in  numerical  ocean  models  than  strict  temporal  accuracy.  The  implicit  treatment 
of  the  vertical  diffusion  terms  should  be  reasonably  accurate  when  the  diffusivities  are  small  and 
the  timescale  of  the  mixing  is  long,  and  when  the  diffusivities  become  large,  it  could  be  argued 
that  the  precise  values  of  the  diffusive  fluxes  are  not  very  well  known  anyway. 

The  version  of  SCRUM  we  received  used  the  Crank-Nicolson  temporal  numerical  scheme  for 
vertical  diffusion.  With  the  Crank-Nicolson  scheme,  the  vertical  diffusion  terms  are  calculated  with 
the  values  of  the  field  being  diffused  averaged  between  the  (n)  and  ( n  + 1)  time  levels,  i.e., 


The  Crank-Nicolson  scheme  is  centered  in  time,  which  would  appear  to  be  an  advantage  over  the 
fully  implicit  scheme  that  is  not  centered.  However,  although  the  Crank-Nicolson  diffusion  scheme 
is,  like  the  fully  implicit  scheme,  unconditionally  numerically  stable,  it  has  the  undesirable  property 
that  the  gradient  of  the  field  being  advected  can  be  reversed  if  the  eddy  coefficients  K  are  sufficiently 
large.  This  gradient  reversal  is  caused  by  “overshoot”  of  the  diffusive  fluxes  that  can  occur  if  the 
vertical  diffusion  number  does  not  satisfy  the  constraint 


2  A  tK 


A  z* 


<  1.0 , 
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where  A z  is  the  vertical  grid  spacing.  With  realistic  surface  atmospheric  forcing,  the  K  values  can 
become  very  large  and  (51)  can  be  easily  violated.  The  resulting  vertical  mixing  tends  to  be  noisy 
and  inaccurate.  This  is  the  reason  most  models  use  fully  implicit  vertical  mixing.  With  fully  implicit 
mixing,  the  diffusive  flux  is  calculated  only  from  the  new  values  of  the  field  being  diffused,  and 
diffusive  overshoots  cannot  occur. 

Vertical  mixing  tests  conducted  with  SCRUM  did  indeed  produce  noisy  and  erratic  results. 
This  behavior  was  eliminated  as  the  timestep  was  reduced  to  the  point  where  the  restriction  defined 
by  (51)  was  approximately  obeyed.  We  note  that  the  newer  version  of  SCRUM,  Version  3.0,  has 
been  converted  to  use  fully  implicit  vertical  mixing,  as  is  used  by  the  other  models. 


4.1  Surface  Mixing 

Some  simple  tests  of  surface-mixed-layer  deepening  were  conducted  to  determine  the  amount 
of  surface  mixing  provided  by  the  models.  The  tests  follow  the  idealized  mixing  experiments 
conducted  by  Martin  (1985,  1986).  One  set  of  tests  consists  of  the  erosion  of  stable  thermal 
stratification  by  a  constant  wind  stress  with  zero  surface  buoyancy  flux.  Another  set  of  tests  con¬ 
sists  of  the  shallowing  of  an  initially  deep  mixed  layer  due  to  steady  surface  heating  and  a  constant 
wind  stress.  A  third  set  of  tests  looks  at  convective  deepening  due  to  surface  cooling. 

Since  the  surface  mixing  tests  consider  only  local  mixing,  the  setup  of  the  model  region  is  not 
especially  critical  as  long  as  non-local  effects  (i.e.,  horizontal  variations)  in  the  central  part  of  the 
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domain  remain  small  during  the  period  of  model  integration.  However,  since  the  models  have  been 
configured  to  allow  the  use  of  periodic  boundaries  in  both  the  x  and  y  directions,  the  models  can 
be  set  up  to  behave  exactly  as  a  local,  one-dimensional  (1-D)  model  in  which  the  horizontal 
variations  are  identically  zero. 

For  the  surface  mixing  tests,  the  models  were  set  up  with  a  4  x  4  grid  of  horizontal  points 
(6x6  for  ECOM-si),  with  periodic  boundaries  in  both  horizontal  directions.  This  setup  gives  a 
2x2  horizontal  grid  of  interior  model  points,  which  is  the  minimum  number  of  horizontal  points 
that  the  models  can  be  run  with  without  significant  code  modification.  With  periodic  lateral  boundaries 
in  both  directions  and  with  horizontally  uniform  initial  conditions  and  surface  forcing,  the  ocean 
models  behave  exactly  like  local,  1-D,  mixed-layer  models. 

A  constant-spaced  vertical  grid  was  used  with  fairly  high  1-m  resolution  to  resolve  the  mixing 
accurately.  The  latitude  of  the  model  region  was  taken  to  be  29.91°  N  where  the  inertial  period  is 
24  h  (i.e.,  /=  2jt/(24  h)).  The  internal  timestep  A tt  was  taken  to  be  600  s.  The  background  viscosity 
and  diffusivity  below  the  mixed  layer  were  set  to  0.1  cm2/s  for  all  the  models. 

For  the  wind-deepening  experiment,  the  initial  temperature  was  defined  as  T(z)  =  24  -  0.05  z, 
where  the  depth  z  is  in  meters  and  the  salinity  was  set  to  a  constant  35  psu.  Three  experiments  were 
conducted  with  constant  wind  stress  values  of  1,  4,  and  16  dynes/cm2.  The  surface-mixed-layer 
depths  (SMLD)  after  1,  2,  and  5  d  are  listed  in  Table  4  and  hourly  values  of  SMLD  are  plotted  in 


Table  4  —  Depth  of  Surface  Mixing  for  Deepening  of  the  Surface 
Mixed  Layer  Due  to  a  Constant  Wind  Stress  and  for  the  Shallowing  of 
the  Mixed  Layer  Due  to  a  Constant  (1  dyne/cm2)  Wind  Stress  and 
a  Constant  Surface  Heat  Flux.  Two  MLDs  are  Listed:  the  First  is 
the  Depth  Where  T  is  0.2°C  Less  Than  the  SST,  and  the  Second  is  the 
Depth  Where  the  Vertical  Diffusivity  Drops  Below  1  cm2/s. 


Wind  Stress  =  1 

Heat  Flux  =  0 

4 
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Mixed-Layer  Depth  (m)  After  24  h 
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Mixed-Layer  Depth  (m)  After  48  h 

POM 
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ECOM 
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4/9 
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21/20 
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Fig.  9  —  Surface  mixed-layer  depth  vs.  time  for  the 
models  for  a  constant  wind  stress  that  is  turned  on  at 
t  =  0.  Three  cases  are  shown  for  wind  stress  values  of 
1,  4,  and  16  dynes/cm2.  The  initial  condition  was  zero 
velocity,  a  constant  salinity  of  35  psu,  and  a  linear 
temperature  stratification  of  5°C/100  m,  with  an  SST 
of  24°C.  The  latitude  is  29.91°  N. 


Fig.  9.  Two  different  SMLDs  are  listed  in  Table  4.  The  first  is  the  depth  at  which  the  temperature 
is  0.2°C  less  than  the  sea  surface  temperature  (SST)  (this  is  also  the  definition  of  the  SMLD  plotted 
in  Fig.  9).  The  second  is  the  depth  at  which  the  vertical  diffusivity  in  the  surface  mixed  layer 
(SML)  decreases  to  less  than  1.0  cm2/s.  For  the  wind-deepening  experiments,  these  two  depths  tend 
to  be  quite  similar. 

For  the  shallowing  experiment,  the  initial  temperature  and  salinity  were  taken  to  be  a  constant 
19°C  and  35  psu,  respectively,  and  the  wind  stress  was  taken  to  be  a  constant  1  dyne/cm2.  Three 
experiments  were  conducted  with  surface  heat  fluxes  of  600,  1200,  and  2400  ly/d.  No  penetration 
of  the  surface  heat  flux  below  the  ocean  surface  (i.e.,  from  solar  radiation)  was  employed.  The 
depths  of  surface  mixing  for  the  shallowing  experiments  after  1,  2,  and  5  d  are  listed  in  Table  4. 
As  for  the  wind-deepening  experiments,  two  SMLDs  are  listed:  the  depth  at  which  the  temperature 
becomes  0.2°C  less  than  the  SST  and  the  depth  at  which  the  vertical  diffusivity  decreases  to  less 
than  1  cm2/s.  In  the  case  of  strong  surface  heating,  there  can  be  quite  a  large  temperature  gradient 
within  the  SML;  hence,  the  depth  at  which  the  temperature  decreases  to  less  than  0.2°C  below  the 
SST  may  significantly  underestimate  the  actual  depth  to  which  surface  mixing  is  occurring. 

The  SMLDs  predicted  by  POM,  ECOM-si,  and  SZM  are  generally  quite  similar  for  both  the 
wind-deepening  and  shallowing  experiments.  POM  and  ECOM-si  are  set  up  with  identical  parameters 
for  the  MYL2.5  turbulence  model;  hence,  the  slight  difference  in  MLD  must  be  due  to  the 
differences  in  the  numerical  schemes.  ECOM-si  gave  slightly  greater  SMLDs  for  the  16  dyne/cm2 
wind-deepening  case  when  the  timestep  was  reduced. 

The  SMLDs  predicted  by  SZM  are  a  little  noisier  than  those  predicted  by  POM  and 
ECOM-si  (Fig.  9).  This  is  likely  due  to  the  fact  that  the  MYL2  turbulence  scheme  does  not  use 
prognostic  equations  to  predict  the  TKE  and  the  turbulent  length  scale  L  The  use  of  prognostic 
equations  tends  to  provide  a  smoother  evolution  of  the  turbulence  fields.  It  is  notable,  however,  that 
the  implementation  of  the  MYL2  mixing  scheme  in  SZM  gives  results  similar  to  the  MYL2.5 
scheme  in  these  tests. 

For  the  tests  of  convective  deepening  of  the  SML  due  to  surface  cooling,  conditions  were  as 
for  the  mixed-layer  shallowing  tests  except  that  the  initial  temperature  was  taken  to  be  a  constant 
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TIME  (d) 


Fig.  10  —  Surface  mixed-layer  depth  vs.  time  for  the 
models  for  convective  deeping  due  to  a  surface  cooling 
of  -1000  ly/d.  The  initial  temperature  was  a  constant 
19°C  from  the  surface  to  50  m,  followed  by  a  linear 
decrease  to  18.5°C  at  100  m.  The  surface  wind  stress 
was  1  dyne/cm2. 


19°C  down  to  50  m,  followed  by  a  small,  linear  decrease  to  18.5°C  at  100  m,  and  the  surface 
heat  flux  was  taken  to  be  -1000  ly/d.  As  for  the  shallowing  tests,  the  surface  wind  stress  was 
1  dyne/cm2.  For  this  test,  the  SMLD  was  defined  as  the  depth  at  which  the  vertical  diffusivity 
decreased  below  1  cm2/s. 

The  SMLDs  from  the  convective  deepening  test  are  plotted  in  Fig.  10.  After  24  h,  the  SMLDs 
predicted  by  POM,  ECOM-si,  and  SZM  were  69,  69,  and  67  m,  respectively,  and  after  48  h,  the 
depths  were  82,  83,  and  81  m.  We  might  expect  that  the  TKE  transport  in  the  MYL2.5  turbulence 
model  used  by  POM  and  ECOM-si  would  generate  larger  SMLDs  in  this  convective  case  than  the 
MYL2  turbulence  model  used  by  SZM,  but  as  found  by  Martin  (1985),  the  difference  is  small. 
Because  of  its  formulation,  the  MYL2  turbulence  model  does  not  provide  any  convective  penetra¬ 
tion  (i.e.,  erosion  of  the  stable  stratification  at  the  base  of  the  turbulent  mixed  layer).  Based  on  the 
results  here  and  previous  tests  (Yamada  and  Mellor  1975;  Martin  1985),  the  MYL2.5  turbulence 
model  provides  only  a  small  amount  of  convective  penetration. 


4.2  Mixing  in  the  Bottom  Boundary  Layer 

Local  mixing  of  the  bottom  boundary  layer  (BBL)  due  to  oscillatory  tidal  currents  was  simulated 
by  setting  up  the  models  with  doubly  periodic  boundary  conditions  as  was  done  for  the  surface 
mixing  cases  in  the  previous  section.  A  surface  pressure  gradient  force  was  imposed  to  mimic 
forcing  by  the  barotropic  M2  tide,  which  has  a  period  of  12.42  h. 

The  surface  pressure  gradient  force  was  applied  to  the  models  through  the  term  that  would 
normally  be  used  to  define  the  local  tidal  potential.  (The  tidal  potential  is  the  effective  force 
generated  by  the  sun  and  the  moon  on  the  oceans  to  drive  the  barotropic  ocean  tides  and  is  the 
equilibrium  position  that  the  ocean’s  surface  would  take  for  a  particular  tidal  mode  if  there  were 
no  other  influences.)  The  term  for  the  tidal  potential  was  added  to  all  the  models  by  subtracting  the 
horizontal  gradient  of  the  tidal  potential  from  the  horizontal  gradient  of  the  surface  elevation  where 
the  latter  appear  in  the  momentum  equations. 
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From  the  horizontally  homogeneous  form  of  the  momentum  equations,  it  can  be  shown  that 
the  tidal  currents  generated  by  a  forcing  term  equivalent  to  a  sinusoidal  surface  elevation  gradient 
in  the  x  direction  of  the  form 


dx 


=Ap  sin  (cor) , 


(52) 


where  t,p  is  the  tidal  potential,  Ap  is  the  amplitude  of  the  surface  elevation  gradient,  co  is  the  tidal 
frequency,  and  t  is  the  time,  will  generate  barotropic  currents 

u-Aucos  (cor)  (53) 
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where  the  current  amplitudes,  Au  and  Av,  are 
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This  is  the  classic  problem  of  the  periodic  forcing  of  a  harmonic  oscillator.  It  can  be  seen  that,  for 
this  idealized  problem,  the  amplitude  of  the  current  goes  to  infinity  as  the  frequency  of  the  forcing 
co  approaches  the  local  inertial  frequency/.  For  co  =  1.405  x  10-4  s-1  (the  frequency  of  the  M2  tide), 
/=  2jt/(24  h),  and  Ap  =  10-5,  the  amplitude  of  the  currents  are  Au  =  95  cm/s  and  Av  =  49  cm/s.  Note 
that  this  analysis  has  not  included  bottom  drag,  which  will  tend  to  modify  the  currents  to  some 
degree.  The  value  of  the  bottom  drag  coefficient  (calculated  from  (21)  with  c^mjn  =0.0025  and 
z0-  0.3  cm)  is  0.0061. 

Besides  the  amplitude  of  the  tidal  forcing,  the  depth  of  mixing  will  depend  on  the  Coriolis 
parameter,  the  initial  stratification,  and  the  parameterization  of  bottom  drag.  The  temperature  and 
salinity  were  initialized  as  for  the  wind-deepening  experiments  in  Sec.  4.1,  i.e.,  the  temperature  was 
defined  by  an  SST  of  24°C  with  a  linear  gradient  of  0.05°C/m,  and  5  was  taken  to  be  a  constant 
35  psu.  The  total  depth  was  taken  to  be  100  m,  the  vertical  grid  spacing  was  a  uniform  1  m,  /  was 
set  to  2jt/(24  h),  and  At  was  600  s.  Note  that  the  total  depth  should  not  affect  the  depth  of  mixing 
in  the  BBL  unless  the  mixing  gets  near  the  surface. 

Figure  11  shows  BBL  depth  (BBLD)  versus  time  for  the  models,  calculated  for  three  different 
values  of  the  tidal  forcing  amplitude  Ap:  2.5  x  10 ~6,  5  x  10~6,  and  1  x  10-5.  The  values  of  the  tidal 
current  amplitude  Au  from  (55)  for  these  three  cases  are  24,  48,  and  95  cm/s. 

The  BBLD  in  Fig.  11  is  defined  as  the  depth  at  which  the  temperature  becomes  0.2°C  greater 
than  the  bottom-layer  temperature.  Note  that  with  this  definition,  the  BBLD  for  the  initial  stratification 
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Fig.  11 — Bottom  mixed-layer  depth  vs.  time  for  the 
models  for  an  oscillating  tidal  current  with  a  period  of 
12.42  h.  The  tidal  current  is  driven  by  an  imposed  surface 
pressure  gradient.  Three  cases  are  shown  with  maximum 
amplitudes  of  the  tidal  current  of  24,  48,  and  95  cm/s. 
The  initial  condition  was  zero  current,  a  constant  salinity 
of  35  psu,  and  a  temperature  profile  defined  by  a  uniform 
stratification  of  5°C/100  m  and  a  bottom  temperature 
of  19°C.  The  latitude  is  29.91°  N,  where  /=  2jt/(24  h). 


that  was  used  is  4  m,  as  can  be  seen  in  Fig.  11.  The  BBLD  defined  by  the  temperature  change  of 
0.2°C  was  generally  within  2-3  m  of  the  BBLD  defined  by  the  depth  of  significant  turbulent  mixing 
(defined  as  the  depth  where  K  >  1  cm2/s). 

The  near-surface  currents  produced  by  the  models  in  these  simulations  agreed  very  closely 
with  the  currents  that  were  expected  based  on  (53—56).  Near  the  bottom,  the  currents  are  reduced 
by  the  bottom  drag  condition  imposed  by  the  quadratic  drag  law  used  in  the  models.  It  is,  of  course, 
the  vertical  current  shear  generated  by  the  bottom  drag  that  generates  the  mixing  in  the  BBL  in  this 
problem.  The  turbulent  mixing  starts  at  the  bottom  and  the  depth  of  the  BBL  increases  as  the 
turbulence  erodes  the  density  stratification  from  below. 

The  BBLDs  predicted  by  the  models  are  quite  similar  and  evolve  similarly  in  time.  The 
BBLDs  predicted  by  SZM  increase  relative  to  the  BBLDs  predicted  by  POM  and  ECOM-si  as 
the  magnitude  of  the  forcing  is  increased. 

In  the  initial  simulations  of  tidal  mixing,  a  problem  was  encountered  with  the  BBLDs  predicted 
by  POM.  The  BBLD  was  much  larger  than  the  values  obtained  with  the  other  models,  and  the 
profiles  of  KH  were  very  noisy.  Investigation  of  the  problem  showed  that  it  was  due  to  round-off 
error  in  calculating  the  density  when  using  32-bit  floating  point  arithmetic.  A  comment  within 
subroutine  DENS  of  POM  recommends  that  several  variables  be  made  double  precision  when  using 
32-bit  arithmetic.  When  the  recommended  variables  in  DENS  were  converted  to  double  precision, 
the  mixing  problems  were  eliminated. 

Note  that  the  density  calculation  used  in  POM  differs  from  that  used  in  the  other  models  in 
that  POM  includes  the  effect  of  pressure  on  the  density.  This  aspect  of  the  density  calculation 
appears  to  be  especially  sensitive  to  round-off  error.  If  the  model  density  includes  the  effect  of 
pressure,  this  must  be  accounted  for  when  calculating  the  buoyancy  term  in  the  TKE  equation, 
which  POM  does.  Otherwise,  vertical  turbulent  mixing  will  be  significantly  (and  incorrectly)  retarded 
by  the  vertical  gradient  of  density  due  to  the  pressure. 
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4.3  Formation  of  a  Tidal  Mixing  Front  with  Stratification  Generated  by  Surface  Heating 


The  formation  of  a  tidal  mixing  front  along  an  infinite  coast  was  simulated  by  setting  up  the 
models  in  what  we  call  “pseudo  2-D”  form.  In  this  form,  one  of  the  horizontal  directions  (in  this 
case,  the  y  direction)  is  set  up  with  the  minimum  number  of  interior  gridpoints  (two)  and  with 
periodic  boundaries.  With  this  setup,  the  models  behave  as  if  they  were  2-D  without  the  effort  to 
actually  convert  the  models  to  a  true  2-D  form  (although,  of  course,  such  a  pseudo  2-D  model  is 
not  as  computationally  efficient  as  a  true  2-D  model  would  be).  Note  that  while  the  pseudo  2-D 
setup  does  not  allow  any  variability  in  the  y  direction,  currents  in  the  y  direction  can  develop. 

The  bathymetry  in  the  offshore  direction  consisted  of  a  linear  slope  of  0.000571  from  a  depth 
of  20  m  at  the  coast  to  a  depth  of  200  m  at  a  distance  of  315  km  from  the  coast,  followed  by  a 
constant  depth  of  200  m  out  to  the  open  boundary  at  a  distance  of  406  km  from  the  coast  (Fig.  12). 
The  horizontal  grid  spacing  was  7  km.  For  the  main  set  of  simulations,  20  sigma  layers  were  used 
in  the  vertical  with  a  spacing  in  the  deep  water  of  5  m  at  the  surface  and  with  a  uniform  stretching 
of  the  grid  to  the  bottom.  Simulations  were  also  conducted  with  SZM  with  z-level  grids  of  20  and 
50  levels. 

Rather  than  impose  an  initial  stratification,  the  thermal  stratification  was  allowed  to  develop 
from  an  initially  homogeneous  condition  under  the  influence  of  surface  heating  and  mixing.  The 
solar  flux  ( Qr )  and  the  downward  surface  heat  flux  (Qf,  +  Qe  +  Qs)  were  taken  to  be  450  and  -300 
ly/d,  respectively.  These  fluxes  give  a  net  surface  heating  of  150  ly/d.  A  two-band  approximation 
of  Jerlov  Type  IA  seawater  (Jerlov  1968)  was  used  to  describe  the  solar  extinction 

y(z)  =  0.38ez/(20  m>  +  0.62ez/(°35  m) .  (57) 

This  parameterization  gives  more  solar  penetration  than  usually  occurs  in  coastal  water,  which 
tends  to  be  more  turbid  than  Type  IA.  However,  this  extinction,  together  with  the  surface  heat 
fluxes,  maintains  an  SMLD  of  about  12  m  (due  to  convection  driven  by  the  solar  penetration) 
without  the  need  to  employ  a  surface  wind  stress  to  provide  mixing  of  the  surface  layer. 

The  M2  tide  was  forced  at  the  open  boundary  by  specifying  a  sinusoidally  varying  surface 
elevation  with  an  amplitude  of  1  m  and  a  period  of  12.42  h.  Table  5  lists  results  from  the  simulations. 
The  values  listed  in  Table  5  include  the  internal  mode  (baroclinic)  timestep  At,  the  timestep  for 
the  free  surface  mode  A te,  the  amplitude  of  the  surface  elevation  at  the  coast,  the  maximum  tidal 
current  that  occurs  within  the  model,  and  the  depth  of  the  water  at  the  location  of  the  front.  The 
tidal  elevation  amplitude  predicted  by  the  models  at  the  coast  is  about  3.5  m  (the  bathymetry  is 
such  that  the  tide  amplifies  significantly  over  the  shelf),  and  the  maximum  tidal  currents  are  about 
70  cm/s. 


Fig.  12  —  Schematic  showing  the  model  domain  used 
for  the  problem  of  the  formation  of  a  tidal  mixing  front 
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Table  5  —  Comparison  of  Model  Results  for  the  Formation  of  a  Tidal  Mixing  Front. 
Tidal  Amplitude  at  Open  Boundary  is  100  cm. 


MAX  AMP 

MAX 

DEPTH  AT 

VERTICAL 

NUMBER 

At 

Afe 

AT  COAST 
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73 

96 

mm 
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Figure  13  shows  temperature  sections  at  30  d  from  the  model  simulations  conducted  with 
vertical  grids  of  20  sigma  layers  and  with  At  =  900  s.  The  results  of  the  models  are  similar.  The 
surface  heating  builds  up  a  vertical  temperature  differential  of  about  1.3°C  in  the  offshore  water. 


Fig.  13  —  Temperature  x-z  sections 
at  30  d  for  the  formation  of  a  tidal 
mixing  front  for  (a)  POM,  (b)  ECOM- 
si,  and  (c)  SZM.  The  tick  marks  along 
the  horizontal  axis  denote  the 
gridpoint  locations.  20  sigma  layers 
were  used  in  the  vertical,  with  an 
upper  layer  thickness  of  5  m  in  deep 
water  and  a  uniform  stretching  to  the 
bottom. 
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Mixing  near  the  bottom  due  to  the  tidal  currents  maintains  a  deep  BBL.  A  tidal  mixing  front  forms 
at  a  depth  of  about  90  m  where  the  BBL  merges  with  the  SML.  The  SST  differential  across  the 
front  is  about  0.7°C.  The  front  is  qualitatively  similar  to  the  observation  of  a  tidal  mixing  front  in 
the  Irish  Sea  shown  by  Simpson  and  James  (1986).  Offshore  of  the  front,  the  water  is  stratified,  and 
inshore  of  the  front  the  water  is  vertically  mixed  from  the  surface  to  the  bottom.  There  is  a  cooling 
of  the  SST  as  the  tidal  mixing  front  is  approached  from  its  offshore  side.  However,  there  is  a 
warming  of  the  water  toward  the  coast  on  the  inshore  side  of  the  front  that  is  caused  by  the  fixed 
surface  heat  flux  being  mixed  into  a  water  column  of  decreasing  depth. 

Table  5  lists  model  results  with  baroclinic  timesteps  of  both  900  and  200  s  to  give  some  idea 
of  the  sensitivity  of  the  simulations  to  the  timestep.  Note  that  a  timestep  of  200  s  is  much  smaller 
than  would  normally  be  needed  to  maintain  numerical  stability  on  a  7-km  grid.  ECOM-si  showed 
the  largest  change  due  to  the  reduction  of  the  timestep,  a  10%  increase  in  the  tidal  elevation  and 
currents.  This  increase  may  be  due  to  the  damping  of  surface  waves  by  ECOM-si’s  fully  implicit 
treatment  of  the  surface  waves  at  larger  timesteps,  which  is  discussed  in  Sec.  5.0.  SZM  showed  a 
decrease  in  tidal  amplitude  with  the  smaller  timestep  and  POM  showed  a  small  increase. 

Figure  14  shows  a  comparison  of  results  from  SZM  for  the  tidal  mixing  front  problem 
conducted  with  different  vertical  grids:  (a)  a  sigma  coordinate  grid  with  20  layers  (this  is  the  same 
result  with  SZM  that  appears  in  Fig.  13c),  (b)  a  20-layer  primarily  z-level  grid  with  18  z-levels  and 
two  sigma  layers  at  the  surface  to  take  up  surface  elevation  changes  (this  grid  provides  the  same 


X-AXIS 


Fig.  14  —  Temperature  x-z  sections 
at  30  d  for  the  formation  of  a  tidal 
mixing  front  using  SZM  with 
different  vertical  grids:  (a)  20  sigma 
layers  as  in  Fig.  13,  (b)  20  layers 
with  z-levels,  except  for  two  sigma 
layers  at  the  surface,  and  (c)  50  layers 
with  z-levels,  except  for  two  sigma 
layers  at  the  surface.  The  20-layer 
sigma  and  z-level  grids  in  (a)  and 
(b)  have  the  same  vertical  spacing 
in  deep  water.  The  50-layer  z-level 
grid  in  (c)  has  a  uniform  4-m  vertical 
spacing  and  exactly  resolves  the  shelf 
slope. 
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resolution  as  the  20-layer  sigma  coordinate  grid  in  the  deepest  water),  and  (c)  a  uniform 
(unstretched)  primarily  z-level  vertical  grid  with  50  total  layers  with  the  two  top  layers  being  sigma 
layers. 

The  results  of  the  simulations  with  the  different  vertical  grids  in  Fig.  14  are  generally  similar. 
However,  one  noticeable  difference  between  the  results  is  the  wavy  disturbance  of  the  isotherms  in 
Fig.  14b,  the  simulation  with  the  grid  of  two  sigma  layers  and  18  z-levels.  In  this  simulation,  the 
fitting  of  the  bathymetry  to  the  z-levels  results  in  a  step-like  bottom,  with  several  horizontal 
gridpoints  along  each  step.  A  result  of  the  steps  is  that  the  convergence  of  the  onshore,  barotropic 
tidal  flow  is  concentrated  at  the  faces  of  the  steps.  Hence,  instead  of  the  vertical  flow  varying 
uniformly  in  the  horizontal,  there  is  a  vertical  “jet”  at  the  face  of  each  of  the  steps  that,  since  the 
flow  is  approximately  barotropic,  reaches  all  the  way  up  through  the  water  column.  A  schematic 
that  illustrates  this  is  shown  in  Fig.  15.  As  can  be  seen  in  Fig.  14b,  the  vertical  jets  displace  the 
isotherms  above  each  of  the  steps. 

Better  resolution  of  the  bathymetry,  which  can  be  obtained  by  increasing  the  number  of 
vertical  levels,  will  reduce  this  effect.  In  Fig.  14c,  the  (primarily)  z-level  grid  with  uniform  (4-m) 
spacing  exactly  resolves  the  bottom  slope.  Hence,  the  horizontal  variation  of  the  convergence  of  the 
onshore/offshore  tidal  flow  and  the  vertical  displacement  of  the  isotherms  occurs  smoothly.  Some 
z-level  models  employ  truncation  of  the  bottom  grid  cells  to  the  actual  bathymetry  to  better  resolve 
the  bathymetry  (Wolff  et  al.  1996). 

The  two  simulations  with  SZM  that  were  conducted  with  z-level  grids  resulted  in  slightly  (3%) 
stronger  tides  than  the  simulation  with  sigma  coordinates  (Table  5). 


4.4  Checkerboard  Mixing 

A  phenomenon  that  occurs  with  all  the  models  in  certain  situations  is  differential  vertical 
mixing  at  alternate  horizontal  gridpoints,  a  phenomenon  that  will  be  referred  to  here  as  “checker¬ 
board  mixing.”  This  was  observed  to  occur  in  the  SML  under  conditions  of  surface  heating  and 
light  winds. 


Fig.  15  —  Schematic  showing  an  x-z  section  for  the 
horizontal  convergence  of  a  barotropic  flow  toward 
the  coast  for  a  model  with  a  z-level  vertical  grid  in  which  the 
bathymetry  is  rounded  to  the  nearest  z-level.  The  horizontal 
convergence  is  localized  at  the  front  of  a  step  and  generates 
a  vertical  jet  at  the  front  of  the  step. 
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An  idealized  experiment  will  be  conducted  here  to  illustrate  how  checkerboard  mixing  occurs. 
The  experiment  was  conducted  on  a  horizontal  grid  of  11  x  11  points  with  doubly  periodic  bound¬ 
aries  and  with  20  layers  in  the  vertical.  The  horizontal  grid  spacing  was  5  km,  the  vertical  grid 
spacing  was  5  m,  and  the  latitude  was  taken  to  be  29.91°  N.  The  initial  conditions  were  T  =  20°C 
and  S  =  35  psu  with  the  ocean  at  rest.  The  forcing  was  a  uniform  wind  stress  of  0.5  dynes/cm2  and 
a  uniform  surface  heating  of  150  ly/d  (73  watts/m2).  The  results  presented  here  were  obtained  with 
POM;  however,  the  other  models  show  the  same  behavior. 

Figure  16  shows  contour  plots  of  the  SST  and  SMLD  after  2  d.  The  tick  marks  along  the  axes 
show  the  gridpoint  locations.  What  has  happened  is  that  the  vertical  mixing  has  mixed  deeper  at 
alternate  gridpoints.  The  SST  shows  the  checkerboarding  pattern  since  the  SST  is  warmer  where  the 
SMLD  is  shallower.  The  checkerboarding  occurs  because  of  the  way  vertical  mixing  is  calculated 
by  the  models  on  a  C  grid.  The  staggered  velocities  (Fig.  17)  are  averaged  to  the  grid-cell  centers 
to  calculate  the  vertical  mixing  coefficients.  The  vertical  mixing  of  T  and  S  occurs  at  the  grid-cell 
centers  where  T  and  5  are  defined.  However,  the  vertical  mixing  coefficients  for  momentum  are 
averaged  back  to  the  staggered  velocity  points  to  calculate  vertical  mixing  of  u  and  v. 

The  checkerboarding  pattern  gets  set  up  because  slightly  weaker  mixing  at  a  point  relative  to 
that  at  an  adjoining  point  is  reinforced  on  the  next  timestep.  The  weaker  vertical  mixing  of  T  and 
S  results  in  stronger  stratification  at  that  location,  which  inhibits  mixing  on  the  next  timestep.  The 
vertical  shear  of  the  horizontal  velocity  remains  fairly  uniform  horizontally  because  of  the  horizontal 
averaging  of  the  vertical  eddy  coefficients  for  momentum.  Hence,  the  vertical  stability  is  governed 
by  the  vertical  mixing  of  T  and  S. 

Figure  17  shows  a  schematic  illustrating  the  checkerboarding  situation.  The  vertical  eddy 
coefficients  at  the  grid-cell  centers  alternate  in  magnitude,  and  the  strength  of  the  vertical  stratification 
alternates  accordingly.  When  the  vertical  eddy  coefficients  are  averaged  to  the  staggered  u  velocity 


Fig.  16  —  (a)  SST  and  (b)  SMLD  at  2  d  for  the  simulation 
at  which  the  temperature  is  0.2°C  less  than  the  SST.  The  c 


1  11 


AXIS 

:  checkerboard  mixing.  The  SMLD  is  defined  as  the  depth 
tour  interval  is  0.1°C  for  the  SST  and  2  m  for  the  SMLD. 
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Fig.  17  —  Schematic  to  help  explain  the  cause  of  checkerboard  mixing.  The 
diagram  shows  an  x-z  section  through  the  center  of  a  row  of  grid  cells. 


points  at  the  edges  of  the  grid  cells,  the  horizontal  variations  are  averaged  out.  Hence,  the  vertical 
eddy  coefficients  for  momentum,  and  as  a  result,  the  vertical  velocity  shears  are  fairly  uniform 
in  the  horizontal.  Once  the  checkerboarding  pattern  is  initiated,  it  reinforces  itself  and  becomes 
locked  in. 

If  the  conditions  in  the  idealized  problem  discussed  above  are  kept  perfectly  uniform  in  the 
horizontal,  the  response  remains  horizontally  uniform  with  identical  mixing  at  each  horizontal 
point.  However,  any  horizontal  inhomogeneity  introduced  into  the  problem  can  trigger  the  checkerboard 
mixing.  In  the  case  shown  in  Fig.  16,  the  SST  at  the  center  gridpoint  of  the  domain  was  increased 
by  0.01°C.  With  the  introduction  of  this  inhomogeneity,  the  checkerboarding  initially  spreads 
laterally  in  the  direction  of  the  applied  windstress,  and  within  2  d  covers  the  entire  domain. 

Checkerboard  mixing  patterns  were  observed  in  several  different  simulations  that  were  conducted, 
including  some  preliminary  simulations  of  the  tidal  mixing  front  problem  in  which  a  weak  wind 
stress  was  applied  and  the  solar  penetration  was  turned  off.  However,  discussions  with  other  modelers 
indicate  that  it  has  not  been  widely  noticed.  Whether  or  not  checkerboard  mixing  occurs  depends 
on  a  number  of  factors  including  the  applied  surface  fluxes,  the  degree  of  solar  penetration,  and  the 
grid  spacing.  It  seems  less  likely  to  occur  or  make  itself  known  when  realistic  temporally  and 
spatially  varying  atmospheric  forcing  is  being  used.  Even  when  it  does  occur,  it  might  be  viewed 
primarily  as  an  aesthetic  problem.  The  discussion  of  checkerboard  mixing  was  provided  here  mainly 
to  acquaint  the  reader  with  the  phenomenon  in  the  event  it  should  be  encountered. 


5.0  TEST  SURFACE  WAVE  PROPAGATION 

Wave  propagation  is  an  important  aspect  of  the  dynamics  of  coastal  regions.  The  next  three 
sections  discuss  tests  of  wave  propagation  in  the  models.  These  tests  consider  freely  propagating 
surface  waves,  freely  propagating  internal  waves,  and  coastal-trapped  waves  including  surface  and 
internal  Kelvin  waves  and  barotropic  shelf  waves. 

Note  that,  for  a  particular  application,  the  accuracy  of  propagation  of  a  particular  type  of  wave 
may  or  may  not  be  significant.  However,  an  ability  to  accurately  simulate  the  propagation  of 
surface,  internal,  and  coastal-trapped  waves  is,  in  general,  desirable  for  a  model  that  is  to  be  applied 
to  a  wide  range  of  coastal  processes. 
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For  all  the  wave  tests,  the  waves  are  defined  to  be  propagating  in  the  ^-direction.  Hence,  the 
model  domains  are  periodic  in  x,  with  the  wavenumber  in  the  ^-direction  chosen  so  as  to  fit  an 
integral  number  of  wavelengths. 

The  models  are  initialized  from  the  linear,  analytical  solution  for  the  waves,  and  the  model 
errors  are  calculated  by  comparison  of  the  model  fields  with  the  analytic  solution.  To  optimize 
comparison  of  the  model  simulations  with  the  linear,  analytic  solutions,  the  amplitudes  of  the 
waves  are  made  small  to  minimize  nonlinear  effects,  bottom  drag  is  set  to  zero,  the  diffusion 
coefficients  are  kept  small,  and  the  lateral  boundaries  for  the  coastal  trapped  wave  simulations  are 
made  free  slip. 

Three  aspects  of  wave  propagation  by  the  models  are  investigated:  (1)  phase  speed  error, 
(2)  damping,  and  (3)  distortion  of  the  wave  form, 

The  phase  speed  error  for  the  models  is  calculated  by  finding  the  phase  speed  that  gives  the 
best  fit  of  the  analytic  solution  to  the  model  solution,  i.e.,  the  minimum  root-mean-squared  (RMS) 
difference  or  maximum  correlation  between  the  analytic  and  model  solutions.  Note  that  the  other 
model  errors  are  calculated  with  the  phase  of  the  analytical  solution  matched  to  that  of  the  model 
solution. 

The  damping  is  calculated  by  finding  the  fractional  damping  factor  D  for  the  analytic  solution 
that  gives  the  minimum  RMS  error  between  the  analytic  and  model  fields  (with  the  phase  of  the 
analytical  solution  matched  to  the  phase  of  the  model  solution).  Originally,  the  damping  was 
calculated  by  comparing  the  magnitudes  of  the  maximum  and  minimum  values  of  the  analytic  and 
model  fields;  however,  this  calculation  tended  to  underestimate  the  model  damping  if  there  was 
any  noise  in  the  model  fields.  The  e-folding  damping  timescale  t ^  is  computed  from  the  elapsed 
time  t  and  the  amount  of  damping  as  td  =  -t/ln(D). 

The  distortion  of  the  waves  is  characterized  by  the  correlation  between  the  analytic  and  model 
fields  (again,  with  the  phases  of  the  fields  matched). 

These  three  types  of  error  can  be  considered  to  be  independent.  For  example,  if  the  wave 
propagates  with  a  certain  phase  speed  error,  but  maintains  its  original  amplitude  (corresponding  to 
zero  damping  error),  its  damping  timescale  will  be  infinite.  Similarly,  if  the  wave  propagates  with 
a  certain  phase  speed  error  and  a  certain  amount  of  damping,  but  perfectly  maintains  its  initial 
form,  the  correlation  of  the  model  solution  with  the  analytic  solution  will  be  exactly  one. 

In  addition  to  calculating  the  phase  speed  error,  damping  error,  and  correlation,  the  RMS 
difference  between  the  analytic  and  model  solutions  is  also  reported.  This  RMS  error  is  calculated 
with  the  phase  of  the  analytic  solution  matched  to  that  of  the  model  solution,  but  without  any 
damping  applied  to  the  analytic  field,  i.e.,  the  RMS  error  will  reflect  both  the  damping  and  distortion 
of  the  model  solution,  but  not  the  mean  phase  speed  error. 


5.1  Description  of  Surface  Wave  Tests 

The  propagation  of  plane  parallel  surface  waves  was  tested  using  a  doubly  periodic  (periodic 
in  both  x  and  y )  model  domain.  With  the  waves  propagating  in  the  x  direction,  there  is  no  variation 
of  the  fields  in  y,  i.e.,  the  problem  is  2-D  with  variation  only  in  x  and  z.  Hence,  the  doubly  periodic 
grid  is  used  with  the  minimum  dimension  in  y  (two  interior  gridpoints)  that  can  be  used  with  the 
3-D  model  codes  (this  was  referred  to  earlier  as  a  quasi  2-D  model  grid). 
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The  waves  were  initialized  in  the  models  from  the  analytical  solution  for  linear,  shallow-water 
(i.e.,  with  wavelength  large  relative  to  the  depth),  surface  gravity  waves  propagating  in  the 
^-direction.  The  equations  for  these  waves  are 


du  dt, 

(58) 

<3v 

(59) 

Tr-fu- 

dt  dx  ' 

(60) 

The  analytical  solution  for  the  waves  is 


£  =  A  cos  (kx  - 1 of)  , 

(61) 

CD 

u  -A—-  cos  (kx  —  (x)t)  , 

Kn 

(62) 

f 

V  =Ay—  COS  (kx  -  (Jit)  , 

kH 

(63) 

where  A  is  the  maximum  amplitude  of  the  surface  displacement,  k  is  the  horizontal  wavenumber, 
and  co  is  the  frequency.  Note  that  k  =  2jt/L,  where  L  is  the  wavelength  and  co  =  2jt/P,  where  P  is 
the  wave  period.  The  dispersion  relation  for  these  waves  is 

co  2  =  c20k2+f2, 

(64) 

where  c0  =  ( gH)\  and  the  phase  speed  is 

<0  I  2j2Y 

e*+?J  • 

(65) 

For  wavelengths  less  than  about  100  km,  the  Earth’s  rotation  does  not  have  much  effect  on 
surface  waves  because  the  period  of  the  waves  is  sufficiently  short  relative  to  the  local  inertial 
period  (2rt//)  that  rotational  effects  do  not  have  much  time  to  act.  The  effect  of  the  wave  frequency 
on  the  phase  speed  is  easier  to  see  if  the  expression  for  the  phase  speed  (65)  is  rewritten  as 


c  =  c0 


1 


(66) 


For  the  shorter  waves,  c  =  cQ  and  the  waves  are  approximately  non-dispersive,  i.e.?  the  phase  speed 
is  independent  of  wavelength.  For  long  wavelengths  of  order  1000  km,  the  phase  speed  depends  on 
the  wavelength  and  the  waves  are  dispersive. 
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Tests  were  conducted  both  for  short  waves  ( L  <  100  km)  and  for  long  waves  (. L  >  1000  km) 
that  are  affected  by  rotation.  The  depth  H  of  the  model  region  was  set  to  40.81  m,  which  gives 
c0  -  20  m/s.  T  and  S  were  set  to  constant  values  of  20°C  and  35  psu,  respectively.  (Note  that  the 
particular  (constant)  values  of  T  and  S  that  are  used  do  not  affect  the  propagation  of  surface  waves 
in  the  models.)  Bottom  drag  was  set  to  zero.  Ten  uniformly  spaced  layers  were  used  in  the  vertical. 
Since  the  bottom  drag  is  set  to  zero,  there  are  no  vertical  variations  in  the  flow  and  the  number  of 
vertical  layers  that  is  used  is  arbitrary.  For  the  Asselin-filtered  leapfrog  schemes  of  POM  and  SZM, 
v  was  set  to  0.05,  which  is  the  value  typically  used  in  the  models. 

The  internal  timestep  At  for  the  explicit  surface  wave  propagation  scheme  used  by  POM  was 
(for  most  of  the  tests)  set  to  a  value  commensurate  with  the  need  to  (in  more  general  circumstances) 
maintain  stability  for  the  internal  wave  propagation  and  horizontal  advection  terms  for  a  maximum 
expected  speed  of  about  2  m/s.  This  is  greater  than  the  maximum  speed  of  internal  waves  and  ocean 
currents  in  most  coastal  areas.  The  external  timestep  for  POM  A te  has  to  be  made  sufficiently  small 
to  maintain  stability  for  the  propagation  of  the  surface  waves.  The  restriction  for  POM’s  explicit 
scheme  is  that  A te  must  be  sufficiently  small  that  a  surface  wave  cannot  travel  a  distance  Ax  in  a 
single  timestep. 

The  models  with  an  implicit  treatment  of  the  surface  waves,  ECOM-si  and  SZM,  integrate  all 
their  terms  with  a  single  timestep.  Since  these  models  treat  surface  waves  implicitly,  their  timestep 
is  not  restricted  by  stability  considerations  for  the  surface  waves;  hence,  their  timestep  is  limited 
by  advection  and  internal  wave  propagation,  as  is  the  internal  timestep  of  POM.  Thus,  the  timestep 
for  ECOM-si  and  SZM  was  generally  set  the  same  as  the  internal  timestep  used  for  POM.  Some 
tests  were  conducted  with  ECOM-si  and  SZM  with  smaller  timesteps  to  look  at  the  effect  of 
timestep  on  the  accuracy  of  surface  wave  propagation  for  their  implicit  schemes. 


5.2  Surface  Wave  Propagation  without  Rotation 

The  accuracy  with  which  the  models  propagate  surface  waves  depends  on  the  numerical  scheme 
and  the  temporal  and  spatial  resolution  used.  A  stability  analysis  of  the  numerical  schemes  used  in 
the  models  for  surface  wave  propagation  with/=0  is  presented  in  App.  A.  The  stability  analysis 
predicts  the  damping  and  phase  speed  error  to  be  expected  for  a  particular  numerical  scheme, 
timestep,  grid  resolution,  and  wavelength.  The  errors  predicted  by  the  stability  analysis  agree  well 
with  the  errors  observed  with  the  models. 

The  results  of  the  tests  of  surface  wave  propagation  are  presented  in  Table  6,  which  lists  the 
phase  speed  error,  damping  timescale,  and  the  correlation  and  RMS  error  after  12  h  for  each  of 
the  models  for  different  values  of  grid  resolution,  internal  and  external  timestep,  and  wavelength. 
If  values  of  the  correlation  are  not  listed  in  the  table,  it  is  because  the  waves  became  too  damped 
at  12  h  to  calculate  reliable  values.  Note  that  the  RMS  error  asymptotes  to  0.7071  cm  as  the  model 
solution  becomes  completely  damped. 

The  most  notable  result  in  Table  6  is  how  much  more  accurately  the  explicit  scheme  propagates 
surface  waves,  especially  surface  waves  of  short  wavelength.  The  reason  is  the  much  smaller 
timestep  used  by  the  explicit  scheme  to  resolve  surface  wave  propagation.  The  relatively 
large  timestep  generally  used  with  the  implicit  models  cannot  resolve  the  propagation  of  short 
surface  waves. 


Table  6  —  Comparison  of  Model  Errors  for  Propagation  of  Surface  Gravity  Waves  for/  = 

Correlation  and  RMS  Error  are  at  12  h. 


42 


Martin ,  Peggion,  and  Yip 


CO 

CO 

04 

wo 

o 

o 

o 

o 

o 

o 

CO 

CO 

rH 

Cs 

o 

H 

CO 

CO 

T“H 

03 

00 

O' 

VO 

O' 

O' 

O' 

O' 

O' 

04 

CO 

VO 

O' 

O' 

CO 

5  ^  S 

o 

CO 

VO 

wo 

o 

wo 

o 

o 

o 

o 

CO 

00 

o 

O 

o 

o 

wo 

pi  »:  w 

o 

o 

CO 

VO 

O' 

r- 

O' 

O' 

0- 

o 

CO 

t> 

O'; 

O' 

O' 

CO 

ua 

o 

d 

o 

o 

o 

o 

d 

o 

d 

o 

o 

o 

o 

d 

o 

d 

o 

Cs 

o 

OS 

04 

Os 

os 

OS 

Os 

O' 

VO 

Os 

OS 

pi 

Os 

Os 

Os 

Os 

04 

Os 

00 

00 

Os 

Os 

os 

Os 

OS 

Os 

Os 

Os 

Os 

WO 

Os 

OS 

Os 

03 

O' 

OS 

o 

Os 

os 

OS 

Os 

Os 

OS 

OS 

Os 

Os 

Os 

00 

wo 

OS 

o 

Os 

Os 

Os 

Os 

Os 

Os 

Os 

Os 

os 

Os 

00 

VO 

cs 

d 

d 

d 

o 

o 

o 

o 

o 

d 

d 

d 

d 

o 

00 

Os 

wo 

CO 

rH 

&  s 

o 

o 

o 

wo 

CO 

VO 

O' 

o 

o 

o 

o 

q 

CO 

wo 

04 

OS 

o 

2  5  S' 

o 

r-’ 

t-H 

H 

d 

o 

o 

o 

r-< 

wo 

04 

rH 

rH 

d 

O'* 

<  s  o 

O' 

rH 

rH 

04 

r-H 

Q  H 

Tt 

04 

03 

ua  etf  ^ 

PC  w 


cm  O' 

O  CO  O  00 
O  O  VO  CO 

O  O  6  N  On 


VO  rj-  VO  00 

WO  VO  t-H 
04 


^  \D  CO  't 
W0  VO  rH 

04 


ooooo  ooooo  ooooo 

rH  r-+  T--  r-H  I  OOOOO  OOOOO 

04  04  03  04  04  04  04  04  04  04 


o  o 
o  o 

0404040404  0304030404 


S2SS2  2 

2  2  2  2  2  OOOOO  2  2  2  2  2  O  2 

OOOOO  OUOUU  N  N  N  N  N  On 

Cl,  Or  Ph  Cl.  Cl,  Pj  U  W  m  B3  W  Cfl  M  tfl  M  W  &o 


A  Comparison  of  Several  Coastal  Ocean  Models 


43 


The  accuracy  of  surface  wave  propagation  in  the  implicit  models  can  be  increased  by  reducing 
the  timestep,  as  is  shown  in  the  lower  part  of  Table  6.  However,  because  both  the  external  and 
internal  modes  of  the  implicit  models  are  integrated  with  the  same  timestep,  reducing  the  timestep 
significantly  increases  the  computer  time  required  to  run  these  models. 

For  the  propagation  of  surface  waves  without  rotation,  ECOM-si  and  SZM  have  equivalent 
phase  speed  errors,  but  ECOM-si  has  much  greater  damping  because  of  its  fully  implicit  treatment 
of  the  terms  governing  the  surface  wave  propagation.  Table  6  shows  that  even  with  a  very  small 
timestep  (10  s),  ECOM-si  strongly  damps  short  surface  waves.  The  terms  governing  the  surface 
waves  in  POM  and  SZM  are  centered  in  time,  which  in  itself  does  not  result  in  any  damping.  The 
damping  of  surface  waves  in  POM  and  SZM  is  due  to  the  Asselin  filter  (App.  A). 


5.3  Surface  Wave  Propagation  for  Long  Waves  Affected  by  Rotation 

In  this  section,  the  propagation  of  surface  waves  is  examined  for  f-  0.7292  x  10-4  s-1,  which 
corresponds  to  a  latitude  of  29.91°  N.  The  depth  is  taken  to  be  40.81  m  as  in  the  previous  section. 
The  wavelength  used  in  these  tests  is  1280  km.  Equation  (65)  gives  a  phase  speed  for  these  waves 
of  24.89  m/s,  and  it  can  be  seen  that  rotation  has  increased  the  phase  speed  of  the  waves  over  that 
(20  m/s)  for  the  case  without  rotation.  The  period  of  the  waves  is  14.28  h,  which  is  sufficiently 
long  for  rotational  effects  to  significantly  affect  the  waves. 

Table  7  shows  the  damping  and  phase  speed  error  calculated  for  the  models  for  grid  resolutions 
of  10-40  km  and  timesteps  of  1800-7200  s.  As  was  the  case  with  /=  0,  the  implicit  treatment  of 
the  surface  waves  is  much  more  damping  than  an  explicit  treatment  when  the  timestep  for  the 
implicit  scheme  is  set  by  the  typical  stability  limitations  for  horizontal  advection  and  internal  wave 
propagation. 

The  damping  in  both  POM  and  SZM  is  due  to  the  Asselin  filter.  The  damping  of  POM  is  about 
20  times  smaller  than  that  of  SZM  because  of  the  smaller  timestep  used  for  the  calculation  of  the 
surface  waves  with  the  explicit  treatment  of  the  external  mode. 

The  damping  of  ECOM-si  is  about  10  times  larger  than  that  of  SZM.  The  severe  damping  of 
ECOM-si  in  these  results  is  due  to  the  combination  of  the  fully  implicit  scheme  and  the 
large  timestep.  Note  that  the  results  obtained  for  ECOM-si  in  Table  7  (as  are  all  the  results  for 
ECOM-si  in  this  report)  are  for  the  AB2  treatment  of  the  Coriolis  terms  that  was  implemented  to 
allow  the  use  of  large  timesteps.  The  errors  calculated  with  ECOM-si  were  somewhat  noisy,  and 
the  damping  and  phase  speed  errors  in  Table  7  are  averages  of  values  taken  over  the  first  few 
timesteps.  Calculations  with  the  original  lagged  treatment  of  the  Coriolis  terms  were,  however, 
even  more  noisy  and  sometimes  became  unstable. 

The  phase  speed  error  for  these  long  wavelength  waves  is  less  than  1%  for  POM.  For 
ECOM-si  and  SZM,  the  phase  speed  errors  are  small  for  the  typical  grid  resolutions  used  in  coastal 
modeling,  though  the  errors  become  significant  (1-10%)  at  the  larger  timesteps  that  might  be  used 
if  these  models  were  applied  to  large  domains  with  a  grid  spacing  of  20-40  km. 


5.4  Tidal  Propagation 

The  primary  importance  of  surface  wave  propagation  in  coastal  models  is  usually  to  provide 
an  accurate  simulation  of  tides  and  wind  setup  or  storm  surge.  Since  the  tides  generally  have  fairly 
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long  wavelengths,  the  tidal  propagation  by  the 
models  tends  to  follow  the  behavior  obtained 
in  the  previous  section  for  long  surface  waves 
with  rotation. 


......  . .  ....  In  this  section,  a  sample  tidal  calculation 

Fig.  18  —  Model  domain  for  the  test  problem  of  tides  .  ...  .  ,  , 

in  a  bay  is  used  to  illustrate  the  general  effect  of  the 

temporal  resolution  of  the  models  on  tidal 
prediction  accuracy.  This  calculation  is  pre¬ 
sented  as  a  single  illustration.  A  thorough  study  of  the  accuracy  of  tidal  prediction  with  the  models 
would  involve  investigation  of  a  number  of  variables  including  both  temporal  and  spatial  resolution, 
bathymetry  variability,  and  nonlinear  tidal  generation. 


The  model  domain  for  this  test  problem  is  a  rectangular  bay  240  km  long  and  32  km  wide  with 
a  uniform  depth  of  10  m  (Fig.  18).  The  domain  is  closed  except  for  an  open  boundary  on  the  east 
end.  The  tidal  elevation  is  specified  at  the  open  boundary  (i.e.,  a  “clamped”  elevation  boundary 
condition)  with  a  sinusoidal  variation,  an  amplitude  of  10  cm,  and  a  period  corresponding  to  the 

M2  tide  (12.42  h).  (A  small  amplitude 
is  used  to  reduce  nonlinear  effects.) 


Table  8  —  Comparison  of  the  Models  for  the  Propagation 
of  the  Tide  in  a  Bay.  The  Magnitude  and  Phase  Listed 
in  the  Table  is  that  for  the  Maximum  Elevation  at  the 
West  End  of  the  Bay. 


The  horizontal  grid  resolution  is  4  km 
and  the  vertical  grid  consists  of  10  equally 
spaced  sigma  layers.  T  and  S  are  taken 
to  be  constant.  The  Coriolis  parameter 
is  /=  2:t/(24  h)  and  the  bottom  friction 
coefficient  is  0.0061.  Since  the  magni¬ 
tude  of  the  vertical  mixing  coefficients 
will  affect  the  tidal  calculation  in  these 
models,  the  vertical  mixing  coefficients 
were  set  to  a  constant  100  cm2/s  for  all 
the  models  to  provide  a  more  consistent 
comparison. 

Table  8  lists  the  maximum  tidal 
elevation  and  corresponding  phase  at 
the  west  end  of  the  bay  for  POM, 
ECOM-si,  and  SZM  for  values  of  At 
ranging  from  50  to  3200  s.  The 
maximum  elevation  values  from 
Table  8  are  plotted  in  Fig.  19.  Note  that 
for  POM,  A te  has  been  taken  to  be  equal 
to  At/20. 

For  small  values  of  At,  the  models 
predict  a  similar  magnitude  and  phase 
for  the  elevation  at  the  west  end  of  the 
bay.  For  the  fine  grids  usually  used  in 
coastal  modeling  (Ax  <  1  km),  the  small 
timesteps  that  are  required  (Ar  <  200  s) 
provide  good  resolution  of  the  tidal 
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Fig.  19  —  Maximum  tidal  elevation  at  west  end  of  bay  for 
POM,  ECOM-si,  and  SZM  for  different  values  of  the 
timestep.  The  tide  was  forced  by  specifying  the  elevation 
at  the  east  end  of  the  bay  with  an  amplitude  of  10  cm  with 
a  period  corresponding  to  the  M2  tide. 


period  and  will  generally  produce  fairly  accurate  tidal  predictions,  depending  on  the  accuracy 
requirements  and  the  particular  aspects  of  the  simulation. 

The  simulations  with  POM  maintained  the  best  accuracy  for  the  tidal  amplitude  at  the  larger 
timesteps  (Table  8,  Fig.  19).  SZM  provides  fairly  good  accuracy  for  this  problem  up  to  At  =  800  s. 
For  larger  timesteps,  there  are  significant  increases  in  SZM’s  errors  due  to  poor  resolution  of  the 
tidal  period.  The  tidal  amplitude  predicted  by  ECOM-si  shows  significant  damping  at  the  larger 
timesteps,  which  is  consistent  with  the  results  of  the  other  surface  wave  propagation  tests. 

As  the  grid  spacing  and  timestep  are  increased,  the  implicit  treatment  of  surface  waves  results 
in  increasingly  severe  temporal  truncation  error  for  tidal  simulation.  In  simple  geometries,  the 
temporal  truncation  error  of  the  implicit  schemes  at  large  timesteps  may  dominate  the  tidal  prediction 
error.  However,  in  tidal  simulations  in  realistic  coastal  environments,  spatial  truncation  error  will 
also  become  more  severe  as  grid  resolution  is  decreased  because  of  poor  representation  of  the 
bathymetry  and  the  coastline.  Consideration  of  spatial  truncation  error  tends  to  limit  the  size  of  the 
grid  spacing  used  in  coastal  tidal  simulations. 


6.0  TEST  OF  INTERNAL  WAVE  PROPAGATION 

Internal  wave  propagation  in  the  models  was  tested  by  comparing  the  model  solutions  to 
analytical  solutions  for  small  amplitude  internal  waves  propagating  in  a  flat-bottomed  ocean  with 
a  linear  vertical  density  gradient. 

The  linearized  equations  for  hydrostatic  (long)  internal  waves  propagating  in  the  x-direction 
in  a  salt-stratified  ocean  are 


du  1  dp 

dt  p0  dx  ’ 


(67) 
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r 

(68) 

II 

1 

(69) 

dS  N2 
dt-feW’ 

(70) 

du  dw 

~  +  -z~  =  0  > 
dx  dz 

(71) 

where  p  is  the  density,  p0  is  a  reference  density,  5  is  the  salinity,  N  =  (figdS/dz)1/2  is  the 
Brunt-Viasala  frequency,  and  (3  is  the  coefficient  of  expansion  for  salinity.  (Salt  stratification, 
rather  than  thermal  stratification,  was  used  for  this  test  because  of  the  high  linearity  between 
salinity  and  density  changes  at  constant  temperature,  i.e.,  (3  is  fairly  constant.)  The  key  to  linearizing 
the  internal  wave  equations  is  taking  the  density  stratification  in  the  salinity  equation  to  be  fixed, 
which  requires  that  the  amplitude  of  the  internal  waves  be  relatively  small  so  that  the  basic 
stratification  is  not  significantly  changed  by  the  propagation  of  the  waves. 

The  analytical  solution  to  the  above  equations  for  freely  propagating  internal  waves  for  a  flat 
bottom  ( H  =  constant)  and  a  linear  density  stratification  (N  =  constant)  is 


\xtrn 

u  -j-~  cos  (mz)  cos  (kx  -  atf)  , 

(72) 

v -A~^  cos  (mz)  sin  (kx  -  cot)  , 

(73) 

w  =  Aw  sin  (mz)  sin  (kx  -  cot)  , 

(74) 

N 2 

S'  =A-^-  sin  (mz)  cos  (kx  -  cot)  , 

(75) 

where  A  is  the  maximum  amplitude  of  the  vertical  displacement,  S'  is  the  deviation  of  the  salinity 
from  the  background  stratification,  m  =  mt/H  is  the  vertical  wavenumber,  and  n  is  the  vertical  mode 
number.  The  dispersion  relation  for  the  internal  waves  is 

»2=*2*V. 

mL 

(76) 

Hence,  the  phase  speed  is 

to  (n2  f2)2 

(77) 
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Equations  (76)  and  (77)  indicate  that  rotational  effects  tend  to  increase  the  frequency  and  phase 
speed  of  internal  waves. 

If  we  have 


(78) 


the  Earth’s  rotation  will  not  significantly  affect  the  propagation  of  the  internal  waves  and  the  waves 
will  be  nondispersive.  For  the  long  internal  waves  being  considered  here  (L  about  10  km),  this 
condition  may  or  may  not  be  satisfied,  depending  on  the  particular  values  of  the  parameters  in  (77). 
Slower  phase  speeds  and  longer  wavelengths  both  increase  the  period  of  internal  waves  and  increase 
the  time  during  which  rotational  effects  can  act  to  modify  the  waves.  Equations  (76)  and  (77)  can 
be  combined  to  give 


1 

N 
c  =  — 
m 

which  makes  it  clear  that  the  influence  of  the  Earth’s  rotation  on  the  propagation  of  internal  waves 
depends  on  the  ratio  of  the  inertial  and  internal  wave  frequencies. 

For  the  tests  of  internal  wave  propagation  conducted  with  the  models,  the  depth  H  was  taken 
to  be  40  m,  the  temperature  was  taken  to  be  a  constant  20°C,  and  a  linear  salinity  stratification  was 
used  with  a  surface  salinity  of  30  psu  and  a  bottom  salinity  of  40  psu.  With  this  stratification, 
N  =  0.0426  s-1.  The  background  vertical  viscosity  for  momentum  was  set  to  0.1  cm2/s.  The  back¬ 
ground  vertical  diffusivity  for  salt  was  set  to  a  lower  value  of  0.001  cm2/s  so  as  not  to  significantly 
diffuse  the  ambient  salinity  stratification.  Bottom  drag  was  set  to  zero.  Sixteen  layers  were  used  in 
the  vertical  with  a  uniform  spacing. 

Tests  were  conducted  both  with  /=  0  and  with  /=  2jt/(24  h).  Table  9  lists  the  analytical  phase 
speeds  for  the  internal  waves  for  various  wavelengths  and  vertical  modes  for  the  constant  stratification 
used  in  the  tests.  Phase  speeds  are  given  for  /  =  0  and  =  2n:/(24  h). 

As  a  matter  of  interest,  non-hydrostatic  phase  speeds  are  also  listed  in  Table  9.  These  internal 
wave  phase  speeds  are  calculated  without  making  the  hydrostatic  approximation,  which  is  used  by 
all  the  models  being  tested  here.  Table  9  illustrates  that  the  hydrostatic  phase  speeds  are  the  same 
as  the  non-hydrostatic  phase  speeds  for  the  longer  internal  waves,  and  are  only  slightly  in  error  for 
1-km  wavelength  waves. 


1- 


/ 


to 


2 


(79) 


6.1  Test  Internal  Wave  Propagation  without  Rotation 

Tables  10  and  11  list  the  errors  calculated  for  the  propagation  of  first  and  fourth  mode  internal 
waves  of  various  wavelengths  with  /  =  0.  The  horizontal  resolution  of  the  waves  ranges  from  64  to 
4  points  per  wavelength.  Since  there  are  a  total  of  16  vertical  layers,  the  mode  1  waves  have 
16  points  per  mode  in  the  vertical  and  the  mode  4  waves  have  4  points  per  mode. 

The  errors  for  POM  and  SZM  are  similar.  Since  these  models  use  basically  the  same  numerical 
scheme  for  the  baroclinic  part  of  the  equations  (leapfrog  in  time  and  second-order,  centered  spatial 
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Table  9  —  Analytical  Phase  Speed  for  Hydrostatic  and  Non-Hydrostatic 
Internal  Waves  for  a  Depth  of  40  m  and  a  Linear  Salinity  Stratification 
with  a  Brunt-Vaisala  Frequency  of  0.0426  1/s 


PHASE  SPEED  (cm/s) 

HYDROSTATIC 

WAVES 

NONHYDROSTATIC 

WAVES 

MODE 

WAVELENGTH  (km) 

/=  o 

2jt 

/_24h 

/=  o 

,  2jt 
;“24h 

1 

64 

54.25 

91.82 

54.25 

91.82 

32 

65.69 

65.69 

16 

57.33 

57.32 

8 

55.04 

55.03 

4 

54.45 

54.24 

54.44 

2 

54.30 

54.21 

54.26 

1 

54.26 

54.08 

54.09 

2 

64 

27.13 

78.88 

27.13 

78.88 

32 

45.91 

45.91 

16 

32.84 

32.84 

8 

28.66 

28.66 

4 

27.52 

27.12 

27.52 

2 

27.22 

27.22 

1 

27.15 

27.13 

4 

64 

13.56 

75.30 

13.56 

75.30 

32 

39.44 

39.44 

16 

22.95 

22.95 

8 

16.42 

16.42 

4 

14.33 

13.56 

14.33 

2 

13.76 

13.76 

1 

13.61 

13.61 

differences),  the  errors  would  be  expected  to  be  similar.  The  magnitude  of  the  errors  increases  as 
the  spatial  resolution  of  the  waves  is  reduced. 

The  phase  speed  and  damping  errors  for  ECOM-si  are  similar  to  those  for  POM  and  SZM.  The 
correlation  errors  for  ECOM-si  are  also  similar  at  the  longer  wavelengths,  but  become  larger  at 
the  shorter,  less  well-resolved  wavelengths,  which  indicates  more  distortion  of  the  waveform.  The 
spatial  treatment  of  the  baroclinic  mode  in  ECOM-si  is  the  same  as  that  in  POM  and  SZM.  However, 
the  temporal  treatment  is  different  since  ECOM-si  lags  the  horizontal  pressure  gradient  and  advection 
terms  in  time.  Even  though  the  internal  waves  are  fairly  well  resolved  in  time  by  the  200-s  timestep, 
ECOM-si  appears  to  suffer  some  effects  of  temporal  truncation  error  when  the  internal  waves  are 
less  well-resolved  by  the  grid. 

Figures  20  and  21  show  plots  of  the  initial  salinity  anomaly  S'  and  the  anomaly  from  the 
models  at  48  h  for  the  mode  1  internal  waves  for  the  64-  and  8-km  wavelength  tests.  The  anomalies 
for  POM  and  SZM  in  Fig.  20  for  the  64-km  waves  show  good  symmetry  both  horizontally  and 


Table  10  —  Comparison  of  Model  Errors  for  Propagation  of  First  Mode  Internal  Gravity  Waves  for 
/=  0.  Correlation  and  RMS  Error  are  at  48  h.  The  RMS  Error  Asymptotes  to  0.1768  psu  as  the  Waves 
Become  Fully  Damped. 
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Table  11  —  Comparison  of  Model  Errors  for  Propagation  of  Mode  4  Internal  Gravity  Waves  for /=  0. 
Correlation  and  RMS  Error  are  at  48  h.  The  RMS  Error  Asymptotes  to  0.1768  psu  as  the  Waves  Become 
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Fig.  20  —  Salinity  anomaly  S'  for  propagation  of  internal  waves  of  64-km  wavelength. 
The  uppermost  plot  shows  the  initial  condition,  and  the  other  plots  show  results 
after  48  h  for  POM,  ECOM-si,  and  SZM.  The  internal  wave  amplitude  was  100  cm 
and  the  timestep  was  200  s. 


vertically.  However,  there  is  a  bit  of  distortion  in  the  plot  for  ECOM-si,  which  is  seen  as  a  slight 
asymmetry  between  the  isohalines  in  the  top  and  bottom  halves  of  the  domain.  This  distortion  is 
more  apparent  in  Fig.  21  for  the  8-km  waves  and  is  due  to  a  slight  decrease  in  the  background 
salinity  in  the  upper  half  of  the  domain  and  an  increase  in  the  bottom  half  (the  salinity  anomalies 
are  computed  based  on  the  original  background  salinity  stratification). 
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Fig.  21  —  Salinity  anomaly  S  for  propagation  of  internal  waves  of  8-km  wavelength. 
The  uppermost  plot  shows  the  initial  condition  and  the  other  plots  show  results  after 
48  h  for  POM,  ECOM-si,  and  SZM.  The  internal  wave  amplitude  was  100  cm  and  the 
timestep  was  200  s. 
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The  modification  of  the  background  salinity  by  ECOM-si  during  internal  wave  propagation  is 
the  result  of  a  numerical  truncation  error  caused  by  a  phase  difference  between  the  velocity  and 
salinity  values  in  the  salinity  advection  term.  An  investigation  of  this  error  is  presented  in  App.  B. 
The  magnitude  of  the  truncation  error  is  approximately  proportional  to  the  timestep  and  the  square 
of  the  wave  amplitude.  For  large  timesteps  and  wave  amplitudes,  the  modification  of  the  background 
stratification  by  internal  wave  propagation  in  ECOM-si  can  be  significant. 

Figure  22  shows  a  plot  of  S'  at  48  h  from  experiments  with  ECOM-si  similar  to  those  shown 
in  Figs.  20  and  21,  but  with  the  timestep  reduced  from  200  to  50  s.  These  plots  show  better 
symmetry  and  the  modification  of  the  ambient  stratification  is  reduced.  Table  10  lists  the  errors  for 
this  experiment  that  show  the  reduction  of  the  correlation  and  RMS  errors  brought  about  by  the 
reduced  timestep. 

The  bottom  of  Table  10  also  shows  errors  for  the  propagation  of  mode  1  internal  waves  with 
POM  and  SZM  with  At  reduced  to  50  s.  The  errors  are  not  much  different  from  those  obtained 
with  the  larger  timestep  except  that  the  damping  is  reduced. 


6.2  Test  Internal  Wave  Propagation  with  Rotation 

Table  12  lists  the  errors  computed  for  the  propagation  of  first-mode  internal  waves  with  a 
wavelength  of  32  km  for/=  2jt/(24  h).  The  period  of  these  waves,  13.53  h,  is  sufficiently  long  for 
the  Earth’s  rotation  to  accelerate  motions  normal  to  the  direction  of  the  waves  and  to  modify  the 
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Fig.  22  —  Salinity  anomaly  S'  for  propagation  of  internal  waves  of  64-  and  8-km  wavelengths 
with  ECOM-si  with  a  timestep  of  50  s.  The  internal  wave  amplitude  was  100  cm. 


Table  12  —  Comparison  of  Errors  in  the  Salinity  Anomaly  Field 
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wave  propagation.  The  phase  speed  of  the  waves  is  65.69  cm/s,  about  20%  larger  than  the  phase 
speed  of  54.25  cm/s  for  the  case  with  /=  0.  The  horizontal  resolution  of  the  waves  ranges  from  32 
to  4  points  per  wavelength.  Since  there  are  a  total  of  16  vertical  layers,  the  waves  have  16  points 
per  mode  in  the  vertical. 

The  phase  speed,  damping,  and  RMS  errors  for  the  waves  are  similar.  ECOM-si  has  the  highest 
correlation  errors,  and  slightly  higher  RMS  errors  than  POM  and  SZM.  Some  of  ECOM-si’s  correlation 
and  RMS  error  may  be  due  to  the  modification  of  the  ambient  stratification  caused  by  internal  wave 
propagation,  which  was  discussed  in  the  previous  section. 


7.0  TEST  PROPAGATION  OF  COASTAL-TRAPPED  WAVES 

Coastal-trapped  waves  are  an  important  part  of  the  dynamics  of  coastal  regions.  In  this  section, 
the  accuracy  with  which  the  models  propagate  several  types  of  coastal-trapped  waves  is  examined. 
The  types  of  waves  considered  are  (1)  barotropic  Kelvin  waves,  (2)  baroclinic  Kelvin  waves,  and 
(3)  barotropic  shelf  waves.  The  basic  question  is  how  well  the  models  propagate  these  kinds  of 
waves  with  a  particular  temporal  and  spatial  resolution. 

In  all  of  these  experiments,  the  models  are  configured  in  a  periodic  channel  with  boundaries 
at  y  =  0  and  y  =  Ly.  The  along-channel  wavenumber  k  is  chosen  so  as  to  fit  a  single  wavelength  in 
the  along-channel  direction.  As  with  the  surface  and  internal  waves,  we  look  at  three  aspects  of  the 
wave  propagation:  phase  speed  error,  damping,  and  distortion  of  the  waveform.  The  wave  propagation 
errors  are  calculated  in  the  same  way  as  the  errors  for  the  surface  and  internal  waves. 


7.1  Barotropic  Kelvin  Waves 

The  equations  governing  linear,  barotropic  Kelvin  waves  propagating  in  the  x-direction  along 
a  coast  are 


du  , 

(80) 

9»  .  K 

rr-fu-g^- 

(81) 

(82) 

dt  ( dx  dyj 

The  ocean  depth  H  is  taken  to  be  constant.  A  boundary  condition  required  for  Kelvin  waves  is  that 
the  cross-shore  velocity  v  is  zero  at  the  coast. 

The  solution  of  this  set  of  equations  for  a  barotropic  Kelvin  wave  is 

£  =  Ae~y/r  cos  (he  -  o it)  ,  (83) 

u  =A  e~y,r  cos  (kx  -  0)t)  , 


(84) 
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v  =  0  ,  (85) 

where  A  is  the  maximum  amplitude  of  the  surface  displacement  at  the  coast  and  r  =  c/f.  The 
dispersion  relation  for  the  surface  Kelvin  waves  is 

to  =  c0k ,  (86) 


where  c0  =  (gH)2.  Hence,  the  phase  speed  of  the  surface  Kelvin  waves  is 

co 

c  =  j=c0.  (87) 

Some  notable  aspects  of  the  Kelvin  wave  solutions  are:  (a)  the  cross-shore  velocity  v  is  zero 
everywhere  (in  the  equation  for  v,  the  Coriolis  term  is  exactly  balanced  by  the  cross-shore  pressure 
gradient),  (b)  the  amplitude  of  the  Kelvin  wave  decays  exponentially  away  from  the  coast  with  an 
e-folding  scale  r,  which  is  referred  to  as  the  Rossby  radius  of  deformation,  (c)  the  cross-shore  scale 
of  the  waves  r  is  relatively  large,  i.e.,  for  water  depths  exceeding  30  m,  r  >  100  km,  (d)  the  speed  of 
the  waves  is  that  for  long  surface  gravity  waves  in  the  absence  of  rotation  and  does  not  depend  on 
/,  (e)  the  waves  are  nondispersive,  and  (f)  the  waves  propagate  only  in  one  direction  along  the 
coast,  with  the  coast  on  the  right  in  the  Northern  Hemisphere  and  on  the  left  in  the  Southern 
Hemisphere. 

For  the  tests  of  surface  Kelvin  wave  propagation  in  the  models,  the  latitude  was  taken  to  be 
29.91°  N  and  the  depth  was  taken  to  be  40.81  m.  These  values  give  c  =  20  m/s  and  an  offshore  scale 
for  the  waves  of  r  =  275  km.  The  along-channel  wavelength  of  the  waves  was  taken  to  be  1280  km. 
The  along-channel  wavelength  was  taken  to  be  relatively  large  to  reduce  the  effect  of  spatial 
truncation  error  in  the  along-channel  direction.  The  width  of  the  channel  was  taken  to  be  640  km, 
which  is  about  2.3  Rossby  radii  (the  Kelvin  wave  propagation  by  the  models  is  not  sensitive  to  the 
width  of  the  channel). 

The  results  of  the  tests  of  surface  Kelvin  wave  propagation  by  the  models  are  presented  in 
Table  13.  Three  main  cases  were  run  with  horizontal  grid  resolutions  of  10,  20,  and  40  km  and 
corresponding  values  of  At  of  1800,  3600,  and  7200  s.  Figure  23  shows  the  analytical  solution  for 
the  surface  elevation  and  velocity  vectors  at  48  h  and  Fig.  24  shows  the  surface  elevation  from  the 
models  at  48  h  for  the  case  with  Ax  =  10  km  and  At  =  1800  s. 

The  relative  results  of  the  models  are  similar  to  those  obtained  for  the  propagation  of  surface 
gravity  waves  in  Sec.  5.1.  The  phase  speed  errors  for  POM  and  SZM  are  similar  to  those  in  Table  7 
for  surface  wave  propagation  in  a  rotating  system  (with  the  same  parameter  values).  The  damping 
timescale  of  the  Kelvin  waves  for  POM  and  SZM  is  about  50%  longer  than  that  of  the  corresponding 
freely  propagating  surface  waves  in  Table  7.  POM’s  explicit  treatment  of  surface  waves  propagates 
the  surface  Kelvin  wave  very  accurately  with  little  phase  speed  error  and  only  weak  to  moderate 
damping. 

The  large  timestep  used  with  the  implicit  schemes  of  ECOM-si  and  SZM  results  in  strong 
damping  of  the  Kelvin  waves,  with  ECOM-si ’s  damping  by  its  fully  implicit  treatment  of  the  free 
surface  being  especially  strong.  Figure  24  shows  that  the  amplitude  of  the  Kelvin  wave  propagated 
by  ECOM-si  (with  At  =  1800  s)  is  almost  completely  damped  at  48  h.  The  damping  in  SZM,  as  in 


Table  13  —  Comparison  of  Errors  in  the  Surface 
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Fig.  23  —  Analytical  solution  for  a  surface  Kelvin  wave 
at  48  h:  (a)  surface  elevation  in  centimeters  and 
(b)  horizontal  velocity  vectors.  The  scaling  arrow  in  the 
velocity  vector  plot  indicates  a  speed  of  0.5  cm/s. 


Fig.  24  —  Comparison  of  surface  elevation  (cm)  from 
the  models  at  48  h  for  the  propagation  of  a  surface 
Kelvin  wave:  (a)  POM,  (b)  ECOM-si,  and  (c)  SZM. 
The  horizontal  grid  spacing  was  10  km  and  the  baroclinic 
timestep  was  1800  s. 
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POM,  is  due  to  the  Asselin  filter,  but  is  much  larger  with  SZM  than  POM  with  the  same  filter 
coefficient  (v  =  0.05)  because  of  SZM’s  much  larger  barotropic  timestep. 

The  implicit  schemes  also  have  significant  phase  speed  errors.  We  note  that  SZM’s  phase  speed 
error  in  Table  13  takes  a  large  jump  when  the  timestep  is  increased  from  3600  to  7200  s.  This  is 
due  to  poor  resolution  of  the  inertial  period  with  At  =  7200  s,  which  adds  to  the  gravity-wave 
propagation  error.  Since  SZM  uses  leapfrog  and  ECOM-si  uses  a  two-time-level  temporal  scheme, 
SZM’s  timestep  is  effectively  twice  as  large  as  ECOM-si ’s  for  the  same  value  of  At,  which  means 
that  SZM  feels  the  effect  of  poor  resolution  of  the  inertial  period  sooner  than  ECOM-si  as  At  is 
increased. 

The  bottom  of  Table  13  shows  errors  from  additional  runs  made  with  ECOM-si  and  SZM  for 
the  case  Ax  =  40  km,  with  At  reduced  from  7200  s  to  360  s,  i.e.,  to  the  value  of  A te  used  by  POM. 
With  this  large  reduction  of  the  timestep,  the  errors  are  significantly  reduced,  as  expected.  With 
At  =  360  s,  SZM’s  errors  are  about  the  same  as  POM’s.  With  At  =  360  s,  ECOM-si’s  phase  speed 
error  is  small,  but  the  damping  by  the  fully  implicit  treatment  of  the  surface  waves  is  still  high. 

From  these  results,  it  can  be  summarized  that  error  in  propagating  long-wavelength  surface 
Kelvin  waves  with  these  models  on  coarse  grids  will  be  due  mainly  to  temporal  truncation  error. 


7.2  Internal  Kelvin  Waves 

The  linearized  equations  for  long  (hydrostatic),  internal  Kelvin  waves  propagating  along  a 
coast  in  the  x-direction  in  a  linearly  salt-stratified  sea  (with  temperature  constant)  are 


du  1  dp 

8t  ~fV~  Po  dx  ’ 

(88) 

dv  1  dp 

dt~  P0dy’ 

(89) 

dp 

ir-pS’ 

(90) 

dS  N2 
dt~pgW’ 

(91) 

du  dv  dw 

—  +  —  -f  —  —  0  . 

dx  dy  dz 

(92) 

The  analytical  solution  of  these  equations  is 

u=  A  e~y^r  cos  (m  z)  cos  (Ax  -  cot)  ,  (93) 


v  =  0  , 


(94) 
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w  =  Awe  ylr  sin  ( mz )  sin  ( kx  -  cot) 


N ^  , 

S'  =  A  —  e  yr  sin  (mz)  cos  (kx  -  oot)  , 


(95) 

(96) 


where  A  is  the  maximum  amplitude  of  the  vertical  displacement  and  S'  is  the  deviation  of  the 
salinity  from  the  background  stratification.  The  dispersion  relation  for  the  Kelvin  waves  is 


co  = 


Nk 

m 


(97) 


where  m  =  rm/H  is  the  vertical  wavenumber  and  n  is  the  vertical  mode  number.  Hence,  the  phase 
speed  is 


to  _N_NH 
km  nn 


(98) 


The  properties  of  internal  Kelvin  waves  are  similar  to  the  properties  of  surface  Kelvin  waves: 
(a)  the  cross-shore  component  of  the  velocity  is  zero,  (b)  the  amplitude  of  the  waves  decreases 
exponentially  away  from  the  coastal  boundary,  (c)  the  phase  speed  of  internal  Kelvin  waves  is  the 
same  as  for  freely  propagating  internal  waves  in  the  absence  of  rotation  and  does  not  depend  on 
/,  (d)  the  waves  for  a  particular  vertical  mode  are  nondispersive,  and  (e)  the  waves  propagate  only 
in  one  direction,  with  the  coast  on  the  right  in  the  Northern  Hemisphere  and  on  the  left  in  the 
Southern  Hemisphere. 


The  cross-shore  scale  of  the  internal  Kelvin  waves,  which  derives  from  the  need  to  balance  the 
Coriolis  term  and  the  cross-shore  pressure  gradient  in  the  v-equation,  is 


r  = 


co  c 

Try 


(99) 


For  internal  waves,  the  distance  r  is  referred  to  as  the  internal  Rossby  radius  of  deformation.  Since 
the  phase  speed  of  internal  waves  is  much  smaller  than  the  phase  speed  of  surface  waves,  the  cross¬ 
shore  length  scale  of  the  internal  Kelvin  waves  is  much  smaller  than  that  for  surface  Kelvin  waves. 

For  the  model  tests  of  internal  Kelvin  wave  propagation,  the  parameters  used  were  the  same 
as  those  used  in  the  internal  wave  tests,  i.e.,  the  depth  H  was  40  m,  temperature  was  a  constant 
20°C,  and  a  linear  salinity  stratification  was  imposed  with  a  surface  salinity  of  30  psu  and  a  bottom 
salinity  of  40  psu.  The  latitude  was  29.91°  N.  With  these  values,  N  =  0.0426  1/s,  c  =  54.3  cm/s,  and 
r  =  7.46  km.  The  along-channel  wavelength  of  the  waves  was  taken  to  be  64  km.  The  amplitude  A 
of  the  Kelvin  waves,  which  is  the  maximum  vertical  displacement  of  the  isohalines  at  the  coast,  was 
initialized  to  a  value  of  1  m. 

With  the  above  parameters,  the  phase  speed  of  the  internal  Kelvin  waves  is  54.25  cm/s.  This 
is  the  same  as  the  phase  speed  of  freely  propagating  internal  waves  in  the  absence  of  rotation, 
and  is  20%  slower  than  the  speed  of  freely  propagating  internal  waves  with/=  2/jt/(24  h)  (Table  9). 

Table  14  shows  errors  for  the  propagation  of  internal  Kelvin  waves  with  the  models  for  various 
values  of  horizontal  grid  resolution  and  timestep.  These  errors  can  be  compared  with  those  in 


Table  14  —  Comparison  of  Errors  in  the  Salinity  Anomaly  Field  S'  for  the  Propagation  of  First-Mode 
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Table  12  for  freely  propagating  internal  waves,  which  were  run  with  the  same  parameter  values.  For 
all  the  models,  the  phase  speed  error  and  damping  are  slightly  lower  for  the  Kelvin  waves  than 
for  the  freely  propagating  internal  waves. 

It  must  be  noted  that  the  internal  Kelvin  waves  propagate  quite  well,  even  with  very  coarse 
resolution  of  the  cross-channel  structure.  For  example,  the  case  with  Ax  =  8  km  provides  only  one 
point  within  the  first  Rossby  radius  ( r  =  7.46  km)  of  the  southern  boundary  of  the  channel. 

Figure  25  shows  plots  of  velocity  vectors  and  S'  from  the  analytical  solution  for  the  internal 
Kelvin  waves  at  48  h.  The  plot  of  u  is  for  the  surface  layer  (1.25-m  depth)  and  the  plot  of  S'  is 
for  the  middle  of  layer  8  (18.75-m  depth).  These  plots  illustrate  the  horizontal  and  vertical  structure 
of  a  mode  1  internal  Kelvin  wave. 

Figure  26  shows  S'  from  the  models  at  48  h  for  the  case  with  Ax  =  4  km  and  At  =  800  s.  It  can 
be  seen  that  the  results  from  the  models  are  very  similar  and  agree  well  with  the  analytical  solution 
(Fig.  25c),  except  for  a  slight  phase  lag  and  a  small  amount  of  damping. 

7.3  Barotropic  Shelf  Waves 

Barotropic  shelf  waves  are  also  referred  to  as  topographic  waves  or  topographic  Rossby  waves. 
Their  propagation  depends  on  the  Earth’s  rotation  and  the  change  in  the  absolute  vorticity  of  the 
water  column  as  the  bottom  depth  changes.  Hence,  rotation  and  changes  in  bottom  depth  are  both 
required  for  their  existence.  A  distinction  between  these  waves  and  planetary  Rossby  waves  is  that 
the  existence  of  planetary  Rossby  waves  depends  on  the  latitudinal  variation  of  the  Coriolis  parameter 
rather  than  changes  in  depth. 

The  linearized  equations  for  barotropic  shelf  waves  propagating  in  the  x-direction  along  a 
straight  coast  are  similar  to  those  for  surface  Kelvin  waves,  except  the  depth  H  must  vary  with  y. 
Also,  to  simplify  the  analytic  solution,  the  time  derivative  of  the  surface  elevation  is  scaled  out  of 
the  continuity  equation.  This  is  permissible  because  this  term  is  small  and  is  not  essential  to  the 
physics  of  small-amplitude,  topographic  waves.  The  equations  are  then 


su  at 

Trfv~sTx' 

(100) 

av  at 

(101) 

o=V“-w\ 

dx  dy 

(102) 

The  analytical  solution  of  this  set  of  equations  is  simplified  for  certain  choices  of  the  function 
H(y).  One  possible  choice  is  the  function 


H(y)=H0e2a-y,  (103) 

where  Ha  is  the  depth  at  the  coast  and  a  defines  the  cross-channel  scale  of  the  change  in  water 
depth.  With  the  bathymetry  given  by  (103),  the  analytic  solution  of  (100-102)  is  (LeBlond  and 
Mysak  1978) 
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X-AXIS 


Fig.  25 —  Analytical  solution  for  an  internal  Kelvin  wave  at  48  h:  (a)  horizontal 
section  of  along-channel  velocity  u  at  1.25-m  depth,  (b)  x-z  section  of  u  at  4  km 
from  southern  boundary,  (c)  horizontal  section  of  salinity  anomaly  S'  at  depth  of 
18.75  m,  and  (d)  x-z  section  of  S'  at  4  km  from  southern  boundary.  Contour  interval 
is  1.0  cm/s  for  u  and  0.02  psu  for  S'. 
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Fig.  26  —  Comparison  of  horizontal  section  of  salinity  anomaly  S'  from  models  at  48  h 
for  propagation  of  internal  Kelvin  wave  for  case  with  4-km  grid  spacing  and  800-s 
timestep:  (a)  POM,  (b)  ECOM-si,  and  (c)  SZM.  The  horizontal  section  is  at  the  center 
of  the  eighth  model  layer  at  a  depth  of  18.75  m.  Contour  interval  is  0.02  psu. 
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A  _o 


£  =  —  e  ay  ((ooa  -/&)  sin  (/y)  +  00/  cos  (/y)  cos  (£*  -  atf))  , 
gk 


(104) 

m  =  Ae~ay  (a  sin  (Zy)  +  l  cos  (ly))  cos  (Ax  -  cor) ,  (105) 

v  =Ake~ay  sin  (ly)  sin  (Ax  -  oof) ,  (106) 

where  A  scales  the  amplitude  of  the  waves.  The  dispersion  relation  for  the  shelf  waves  is 

2a  kf 


co  = 


k2  +  l2  +  a2 


(107) 


The  cross-shore  wavenumber  l  is  restricted  to  the  values  /  =  rm/Ly  (where  n  =  1,  2,  3,...  is  the 
cross-channel  mode  number)  so  as  to  satisfy  the  condition  of  zero  flow  normal  to  the  walls  of 
the  channel.  The  phase  speed  is  then 


co  2  af 
*V  +  i2  +  a2' 


(108) 


There  is  a  lot  of  freedom  in  defining  parameters  for  shelf  waves.  The  bathymetry,  channel 
width,  and  along-shore  wavelength  can  all  be  specified  independently,  and  for  each  set  of  these 
values,  there  is  an  infinite  number  of  modes  in  the  cross-channel  direction. 

For  the  model  tests,  the  depth  at  the  coast  H0  was  taken  to  be  20  m  and  the  cross-channel 
e-folding  scale  of  the  bathymetry  was  taken  to  be  (2a)-1  =  40  km.  The  width  of  the  channel  Ly  was 
taken  to  be  64  km.  Hence,  the  depth  in  the  channel  varies  from  20  m  at  the  coastal  boundary  to 

99  m  at  the  offshore  boundary. 

The  models  were  initialized  from  the  analytic  solution  (104-106)  for  a  single  mode  in  the 
cross-channel  direction,  and  for  an  along-channel  wavelength  of  128  km.  These  parameters  give 
a  phase  speed  for  the  shelf  wave  of  36.54  cm/s.  The  amplitude  A  of  the  wave  was  taken  to  be 

100  m2/s,  which  gives  a  maximum  amplitude  of  the  surface  elevation  of  about  0.05  cm  and  maximum 
velocities  of  about  0.5  cm/s.  With  these  small  amplitudes,  the  linear  shelf  wave  can  be  propagated 
quite  accurately  by  the  models. 

Table  15  shows  the  error  in  the  cross-channel  velocity  from  the  models  for  horizontal  grid 
resolutions  ranging  from  2  to  16  km.  These  grid  resolutions  provide  from  64  to  8  points  per 
wavelength  in  the  along-channel  direction  and  from  32  to  4  points  in  the  cross-channel  direction. 
The  surface  elevation  and  velocity  calculated  from  the  analytical  solution  at  48  h  are  shown  in 
Fig.  27,  and  a  comparison  of  the  cross-channel  velocity  fields  from  the  models  at  48  h  for  the  case 
with  Ax  =  8  km  is  shown  in  Fig.  28.  The  cross-channel  velocity  from  the  models  can  be  compared 
with  the  analytically  calculated  cross-channel  velocity  shown  in  Fig.  27. 

The  errors  in  Table  15  increase  as  the  resolution  of  the  wave  decreases,  as  expected.  With 
64  points  per  wavelength,  the  shelf  wave  is  propagated  very  accurately  with  little  phase  error  and 
a  damping  timescale  that  exceeds  2  y.  With  8  points  per  wavelength  and  only  4  points  in  the  cross- 
channel  direction,  the  wave  is  still  propagated  fairly  well,  except  that  the  phase  speed  has  a  large 
error. 
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Fig.  27  —  Analytical  solution  for  a  barotropic  shelf  wave 
at  48  h:  (a)  surface  elevation,  (b)  cross-channel  velocity, 
and  (c)  horizontal  velocity  vectors.  The  contour  interval 
for  the  surface  elevation  is  0.005  cm  and  for  the  cross¬ 
channel  velocity  is  0.05  cm/s.  The  scaling  arrow  in  the 
velocity  vector  plot  indicates  a  speed  of  0.5  cm/s. 


With  ECOM-si,  the  cross-channel  velocity  field  increased  in  amplitude  at  the  beginning  of  the 
shelf  wave  simulations  and  then  decreased  slowly.  The  initial  increase  ranged  from  1.3%  for 
Ax  =  2  km  to  10%  for  Ax  =  16  km.  The  cause  of  the  increase  was  not  determined.  The  amount  of 
the  initial  increase  was  not  affected  by  the  timestep  and  the  initial  fields  used  were  checked  against 
those  used  for  the  other  models.  Note  that  the  damping  timescale  in  Table  15  reflects  only  the  decrease 
in  the  wave  amplitude  after  the  initial  increase.  The  RMS  error  for  ECOM-si  in  Table  15  does,  however, 
reflect  the  initial  increase,  and  this  is  the  reason  it  is  larger  than  the  RMS  error  for  the  other  models. 

Simulation  of  the  shelf  waves  with  ECOM-si  with  the  original,  lagged  treatment  of  the  Coriolis 
terms  was  also  tried.  With  this  treatment  of  the  Coriolis  terms,  the  amplitude  of  the  shelf  waves 
increased  steadily  with  time.  The  increase  at  48  h  was  2.5%  with  Ax  =  2  km  and  Af=400  s  and  was 
9.2%  with  Ax  =  8  km  and  At  =  1600  s. 

All  the  simulations  of  the  barotropic  shelf  waves  discussed  so  far  have  used  sigma  vertical 
coordinates.  Some  additional  simulations  were  conducted  with  SZM  to  investigate  the  propagation 
of  the  barotropic  shelf  wave  using  z-levels.  Figure  29  shows  results  at  48  h  using  SZM  with 
Ax  =  2  km  (64  points  per  wavelength)  and  10  uniformly  spaced  z-levels  (Az  =  10  m).  The  results 
look  quite  noisy  compared  with  the  analytical  solution  in  Fig.  27  and  with  the  sigma  coordinate 
simulations  in  Fig.  28  with  significantly  less  horizontal  and  temporal  resolution.  A  problem  here 
is  that  with  the  bathymetry  truncated  to  the  nearest  model  level  in  SZM,  the  specified  bathymetry 
(103)  is  not  very  well  resolved. 
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X-AXIS 


Fig.  28  —  Comparison  of  the  cross-channel  velocity  fields  at  48  h  from  the  models  for 
the  propagation  of  a  barotropic  shelf  wave:  (a)  POM,  (b)  ECOM-si,  (c)  SZM.  The 
horizontal  grid  spacing  for  the  models  is  8  km,  and  the  barotropic  timestep  is  80  s 
for  POM  and  1600  s  for  ECOM-si  and  SZM.  The  contour  interval  is  0.05  cm/s. 
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X-AXIS 


Fig.  29  —  (a)  Surface  elevation,  (b)  cross-channel  velocity, 
and  (c)  horizontal  velocity  vectors  at  48  h  for  a  shelf 
wave  simulated  by  SZM  with  a  z-level  vertical  grid  with 
Az  =  10  m.  The  horizontal  grid  spacing  is  2  km  and  the 
timestep  is  400  s.  The  contour  interval  is  0.005  cm  for 
the  surface  elevation  and  0.05  cm/s  for  the  cross-channel 
velocity.  The  representative  scaling  arrow  in  the  velocity 
vector  plot  indicates  a  speed  of  0.5  cm/s. 


What  SZM  does  (approximately)  is  propagate  the  shelf  wave  that  is  consistent  with  the 
step-wise  bathymetry  that  is  used  with  the  model  rather  than  the  shelf  wave  derived  for  the  smooth, 
exponentially  varying  bathymetry  that  was  used  to  calculate  the  analytical  solution.  Figure  30 
shows  the  analytical  solution  (see  App  C  for  the  derivation)  for  the  cross-shelf  Velocity  for  the  step¬ 
wise  bathymetry  used  in  the  SZM  model  simulation  in  Fig.  29.  The  analytical  solution  for  the 
step-wise  bathymetry  looks  much  like  the  model  solution  for  the  step-wise  bathymetry. 

Figure  31  shows  results  for  the  cross-channel  velocity  at  48  h  for  simulations  conducted  with 
SZM  with  z-levels  with  finer  vertical  grid  resolution,  i.e.,  with  (a)  30  uniformly  spaced  z-levels 
(Az  =  3.33  m),  (b)  60  uniformly  spaced  z-levels  (Az  =  1.67  m),  and  (c)  51  z-levels  with  a  vertical 
stretching  of  the  grid  employed  so  as  to  exactly  resolve  the  bathymetry  represented  by  (103).  Errors 
for  these  simulations  are  reported  in  Table  16.  From  Fig.  31  and  Table  16,  it  can  be  seen  that  the 
shelf  wave  propagated  on  the  z-level  grid  becomes  more  like  the  analytical  solution  for  the  smoothly 
varying  bathymetry  when  the  true  bathymetry  is  more  accurately  represented. 

The  best  solution  with  the  z-level  grids  is  obtained  with  the  stretched  grid  of  51  points  that  was 
set  up  to  exactly  resolve  the  bathymetry.  However,  the  errors  for  this  case  (Table  16)  are  still  larger 
than  the  errors  obtained  for  the  simulations  with  sigma  coordinates,  and  the  damping  is  signifi¬ 
cantly  greater  than  the  damping  for  the  simulations  with  sigma  coordinates  (Table  15).  A  problem 
with  the  z-level  grid,  besides  difficulty  in  resolving  the  bathymetry,  is  interaction  of  the  horizontal 
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Fig.  30  —  Analytically  calculated  cross-channel  velocity 
field  for  a  barotropic  shelf  wave  propagating  along 
step-wise  bathymetry.  The  step-wise  bathymetry  is  the 
same  as  that  used  by  the  SZM  shelf-wave  simulation 
with  a  z-level  grid  with  Az  =  10  m. 


10  20  30  40  50  60 

X-AXIS 


Fig.  31  —  Cross-channel  velocity  fields  at  48  h  for  the 
propagation  of  a  barotropic  shelf  wave  with  SZM 
for  different  z-level  vertical  grids:  (a)  30  uniformly 
spaced  z-levels,  (b)  60  uniformly  spaced  z-levels,  and 
(c)  51  z-levels  with  variable  spacing  so  as  to  exactly 
resolve  the  exponentially  varying  bathymetry.  The 
horizontal  grid  spacing  for  these  results  is  2  km  and 
the  timestep  is  400  s. 
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Table  16  —  Errors  in  the  Cross-Channel  Velocity  for  the  Propagation 
of  Barotropic  Shelf  Waves  Using  SZM  with  Different  z-Level  Vertical 
Grids.  The  Grid  with  51  z- Levels  is  Stretched  so  as  to  Match  the 
Exponentially  Varying  Bathymetry  Exactly.  Horizontal  Resolution  for 
all  the  Runs  is  64  Points  Per  Wave  with  Ax  =  2  km,  At  =  400  s. 


VERTICAL  GRID 

A  z 
(m) 

SPEED 

ERROR 

(%) 

DAMP 

TIME 

(h) 

COR 

RMS 

ERROR 

(cm/s) 

10  Uniform  z-Levels 

10.00 

1.6 

200 

0.95540 

0.0548 

30  Uniform  z-Levels 

3.33 

-1.0 

300 

0.98777 

0.0339 

60  Uniform  z-Levels 

1.67 

0.5 

330 

0.99862 

0.0240 

51  Stretched  z-Levels 

variable 

-0.5 

330 

0.99998 

0.0225 

currents  with  the  bathymetry  steps  near  the  bottom.  Two  consequences  of  this  interaction,  relative 
to  the  sigma  coordinate  simulations,  are  that  the  damping  is  greater  and  the  horizontal  currents  are 
not  barotropic  near  the  bottom. 


8.0  TEST  FORMATION  OF  UPWELLING/DOWN WELLING  FRONTS 

The  models  were  compared  for  the  formation  of  upwelling  and  downwelling  fronts  using  a 
2-D  model  domain  with  idealized  geometry.  Similar  upwelling  and  downwelling  problems  have 
recently  been  investigated  by  Allen  et  al.  (1995)  and  Allen  and  Newberger  (1996)  with  POM. 

Figure  32  shows  the  model  domain,  which  consists  of  a  symmetric,  2-D  basin  with  gently 
sloping  shelves  on  both  the  east  and  west  sides  out  to  100  km,  steeper  slope  regions  from  100  km 
out  to  200  km,  and  an  interior  of  depth  1000  m.  The  depth  at  the  shelf  break  is  100  m  and  the  total 
width  of  the  basin  is  604  km.  To  avoid  “drying”  of  the  shallowest  points  near  the  coast,  the 
minimum  static  bottom  depth  was  set  to  5  m.  The  gradients  of  the  bottom  slope  along  the  shelf 
and  slope  regions  are  0.001  and  0.009,  respectively. 

The  models  were  run  in  “pseudo  2-D”  mode  for  this  problem,  with  periodic  boundary  conditions 
and  two  interior  points  in  the  along-shore  (y)  direction. 

The  initial  condition  consists  of  a  horizontally  uniform  thermal  stratification  (Fig.  32)  with  the 
ocean  at  rest.  The  initial  SST  is  20°C  and  most  of  the  thermal  stratification  lies  between  50-  and 
100-m  depth.  Salinity  is  set  to  a  constant  35  psu.  The  forcing  consists  of  a  uniform,  along-shore, 
northward-directed  surface  wind  stress  of  2  dynes/cm2.  The  wind  stress  is  ramped  up  linearly  in 
time  over  24  h.  Surface  heat  fluxes  are  zero.  All  the  results  presented  are  at  40  d. 

The  models  were  run  at  two  different  grid  resolutions.  One  set  of  simulations  was  conducted 
with  a  2-km  horizontal  grid  and  50  vertical  layers/levels  with  a  5-m-thick  layer  at  the  surface  in 
the  deep  water  and  uniform  stretching  to  the  bottom.  For  the  second  set  of  simulations,  the  horizontal 
resolution  was  reduced  to  10  km  and  the  vertical  grid  was  reduced  to  20  layers/levels  with  a 
surface  layer  thickness  in  the  deep  water  of  10  m  and  a  uniform  stretching  to  the  bottom. 
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Fig.  32  —  Model  domain  and  initial  temperature  field 
for  upwelling/downwelling  simulations 


Figure  33  shows  temperature  and  along-shore  velocity  at  40  d  for  the  model  simulations  with 
the  2-km  grid  and  50  vertical  sigma  layers.  Table  17  lists  the  minimum  SST  on  the  upwelling 
side,  the  maximum  SST  on  the  downwelling  side,  and  the  maximum  along-shore  currents.  The 
northward  wind  stress  induces  upwelling  on  the  west  side  of  the  domain  and  downwelling  on 
the  east  side.  Strong  along-shore  currents  develop  that  are  in  approximate  geostrophic  balance 
with  the  cross-shore  pressure  gradients.  The  temperature  sections  show  good  general  agreement.  A 
curious  feature  that  appears  in  the  simulations  is  a  sharp  rise  in  the  isotherms  just  seaward  of  the 
shelf  break  on  the  downwelling  side.  This  is  caused  by  vertical  mixing,  which  is  enhanced  at  this 
location  by  the  strong  shears  in  the  along-shore  current. 

The  velocity  sections  show  maximum  along-shore  currents  of  130-140  cm/s  on  the  upwelling 
side  and  190-200  cm/s  on  the  downwelling  side  (Fig.  33,  Table  17).  The  maximum  currents  are 
located  just  seaward  of  the  shelf  break.  As  in  the  experiments  of  Allen  (1995,  1996),  the  along¬ 
shore  current  on  the  upwelling  side  is  a  maximum  at  the  surface  and  decreases  with  depth,  whereas 
on  the  downwelling  side,  the  maximum  current  extends  quite  deep.  The  along-shore  currents  for  the 
models  agree  quite  well. 

Figure  34  shows  temperature  and  along-shore  velocity  for  the  simulations  with  the  10-km  grid. 
The  accuracy  of  the  10-km  simulations  can  be  judged  based  on  the  2-km  results.  The  upwelling 
front  is  more  diffuse  and  the  along-shore  currents  are  weaker  for  the  10-km  simulations  because 
of  the  reduced  grid  resolution  and  the  resulting  stronger  horizontal  mixing.  The  weakest  currents 
in  the  10-km  simulations  are  for  ECOM-si.  The  10-km  simulation  with  ECOM-si  was  rerun  with 
the  timestep  reduced  from  1200  to  200  s  to  see  how  much  difference  might  be  due  to  temporal 
truncation  error.  However,  with  the  smaller  timestep,  the  currents  were  only  slightly  higher  (Table  17). 

A  notable  occurrence  in  all  the  sigma  coordinate  simulations  is  the  warming  of  the  near-surface 
water  above  the  initial  SST  of  20°C  on  the  outer  part  of  the  shelf  on  the  downwelling  side.  This 
warming  cannot  really  be  ascertained  in  Figs.  33  and  34  because  of  the  shortage  of  contour  labels. 
However,  the  warming  can  be  seen  in  Fig.  35a,  which  shows  an  expanded  view  of  the  temperature 
field  in  the  downwelling  area  from  the  10-km  simulation  with  POM.  Also,  the  maximum  SST 
values  at  40  d  for  all  the  runs  are  reported  in  Table  17. 

Since  there  is  no  explicit  heating  in  these  simulations,  we  should  expect  that  the  temperatures 
at  40  d  would  not  exceed  the  initial  maximum  temperature  of  20°C.  However,  the  maximum 
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Table  17  —  Results  from  Upwelling  Downwelling 
Experiments  at  40  d 


UPWELLING 

SIDE 

DOWNWELLING 

SIDE 

Ax 

At 

MIN  SST 

MAX  V 

MAX  SST 

MAX  V 

MODEL 

(km) 

(s) 

(°C) 

(cm/s) 

(°C) 

(cm/s) 

Results  with  Sigma  Coordinate  Grid 

POM 

2 

200 

8.4 

138 

20.5 

202 

ECOM 

2 

200 

8.3 

135 

20.9 

190 

SZM 

2 

200 

7.8 

131 

20.5 

202 

POM 

10 

1200 

8.9 

108 

22.2 

170 

ECOM 

10 

1200 

9.4 

94 

22.4 

143 

SZM 

10 

1200 

8.4 

103 

22.0 

162 

ECOM 

10 

200 

9.7 

96 

22.3 

146 

Results  with  z-Level  Grid 

mm 

2 

1200 

8.4 

mm 

20.0 

200 

m 

10 

1200 

10.2 

mm 

20.2 

156 

temperature  on  the  downwelling  shelf  at  40  d  is  20.5-20. 9°C  for  the  2-km  simulations  and 
22.0-22.4°C  for  the  10-km  simulations  (Table  17).  This  spurious  heating  is  caused  by  the  procedure 
that  was  used  for  horizontal  diffusion  of  T  and  S  in  which  the  horizontal  mean  profile  is  subtracted 
from  the  T  and  S  fields  when  performing  horizontal  diffusion  in  sigma  coordinates  (this  is  discussed 
in  Sec.  2.9).  The  anomaly  from  the  horizontal  mean  temperature  profile  in  the  downwelling  area 
is  shown  in  Fig.  35b  and  is  a  maximum  of  7-8°C  at  the  outer  edge  of  the  shelf.  This  large,  positive 
anomaly  diffuses  heat  in  all  directions  and  the  diffusion  toward  the  shelf  causes  the  warming  of  the 
water  on  the  shelf.  The  spurious  warming  is  reduced  as  the  horizontal  grid  resolution  is  increased 
because  the  horizontal  diffusion  coefficients  are  reduced. 

This  illustrates  a  pitfall  of  subtracting  the  horizontal  mean  profile  from  a  field  when  performing 
horizontal  diffusion  in  sigma  coordinates.  It  might  be  preferable,  instead,  to  subtract  a  smooth 
but  horizontally  varying  field  that  would  be  able  to  account  for  large  departures  from  the  basin 
horizontal  mean.  Such  a  field  would  have  to  be  recomputed  periodically  if  it  changes  significantly 
in  time.  This  procedure  was  not  tried  here,  but  has  been  used  by  other  sigma-coordinate  modelers. 
It  is  noted  that  in  POM,  the  arrays  used  to  hold  the  smoothed  values  of  T  and  S  that  are  subtracted 
from  the  T  and  S  fields  when  performing  horizontal  diffusion  are  most  recently  referred  to  as 
“climate  fields”  (Mellor  1996),  which  suggests  fields  that  have  3-D  structure,  but  are  spatially 
smooth.  A  general  concern  with  such  methods  would  be  to  avoid  spurious  development  of  the  field 
that  might  be  caused  by  “steering”  the  field  via  the  diffusion  term  towards  some  computed  “mean” 
structure,  especially  if  that  mean  is  recomputed  periodically  from  the  evolving  fields. 

The  upwelling/downwelling  simulations  with  2-  and  10-km  horizontal  resolution  were  repeated 
with  SZM  with  a  z-level  vertical  coordinate.  The  z-level  grids  were  set  up  to  have  the  same  vertical 
resolution  in  the  deep  water  as  the  sigma  coordinate  grids.  Figure  36  shows  plots  of  temperature 
and  along-shore  velocity  with  SZM  with  the  z-level  grids,  and  Table  17  lists  the  maximum  and 
minimum  SST  and  the  maximum  along-shore  currents.  With  2-km  horizontal  resolution  and  50  vertical 
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—  Temperature  and  along-channel  velocity  at  40  d  for  the  formation  of  upwelling  and  downwelling  fronts  with  (a)  POM, 
(b)  ECOM-si,  and  (c)  SZM  with  a  horizontal  grid  resolution  of  10  km.  20  sigma  layers  are  used  in  the  vertical. 
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Fig.  35  —  An  expanded  view  of  the  temperature  field  in 
the  downwelling  area  at  40  d  from  the  10-km  simulation 
with  POM.  The  temperature  values  are  shown  in  (a)  and 
the  anomaly  from  the  initial  temperature  profile  is  shown 
in  (b). 


layers,  the  results  with  SZM  with  z-level  and  sigma  coordinates  are  fairly  similar.  The  maximum 
along-shore  currents  are  about  the  same  (Table  17).  The  minimum  SST  on  the  upwelling  side  with 
z-levels  (8.4°C)  is  greater  than  was  obtained  with  SZM  with  sigma  coordinates  (7.8°C),  though  it 
is  similar  to  the  values  obtained  with  POM  and  ECOM-si. 

With  10-km  horizontal  resolution  and  20  vertical  layers,  the  differences  between  the  SZM 
z-level  and  sigma  coordinate  results  are  larger,  e.g.,  the  along-shore  current  and  cooling  on  the 
upwelling  side  show  larger  differences  (Table  17).  Some  of  the  differences  are  due  to  reduced 
vertical  resolution  with  the  z-level  grid  in  shallower  water,  which  provides  reduced  resolution  of 
both  the  model  fields  and  the  bathymetry.  However,  an  experiment  conducted  with  a  10-km  grid 
with  higher  vertical  resolution  showed  that  some  differences  remained,  i.e.,  some  of  the  differences 
are  due  to  other  factors. 

It  can  be  noted  that  the  spurious  warming  of  the  near-surface  temperature  that  occurs  on  the 
downwelling  shelf  with  sigma  coordinates  does  not  occur  when  z-levels  are  used.  The  maximum 
SST  of  20.2°C  that  occurs  in  the  10-km  simulation  with  z-levels  is  probably  due  to  a  small  numerical 
advective  “overshoot”  that  can  occur  at  fronts  with  the  second-order,  centered  advection  schemes 
used  in  these  models  (Sec.  3.2). 

The  trends  seen  in  these  results  illustrate  what  we  might  expect,  that  as  grid  resolution  is 
increased,  differences  between  simulations  conducted  with  sigma  and  z-level  coordinates  tend  to  be 
reduced.  At  coarse  resolution,  the  sigma  and  z-level  coordinates  each  have  particular  advantages 
and  disadvantages. 


9.0  SUMMARY 

Several  ocean  models  that  are  being  used  for  coastal  ocean  simulation  and  prediction,  POM, 
ECOM-si,  SZM,  and  SCRUM  (Version  2.1)  were  tested  for  their  ability  to  simulate  physical  processes 
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of  importance  in  the  coastal  ocean.  The  basic  processes  for  which  the  models  were  tested  included 
advection,  vertical  mixing,  propagation  of  surface,  internal,  and  coastal-trapped  waves,  and  the 
formation  of  upwelling  and  downwelling  fronts. 

Because  of  some  developmental  problems  with  the  version  of  SCRUM  that  we  initially  obtained, 
and  the  later  release  of  an  extensively  modified  SCRUM  code  (Version  3.0)  by  Rutgers,  SCRUM 
was  not  included  in  many  of  the  model  tests  that  were  conducted. 

The  model  tests  revealed  some  particular  problems  with  the  individual  models,  which  will  be 
briefly  summarized  here.  In  many  cases,  these  are  limitations  that  occur  in  certain  situations.  These 
problems  and  limitations  should  be  kept  in  mind  when  applying  the  models  or  when  selecting  a 
model  for  a  particular  application. 

The  forward  time-differencing  scheme  used  by  ECOM-si  suffers  from  significantly  higher 
temporal  truncation  error  than  the  leapfrog  and  Adams-Bashforth  schemes  used  by  the  other  models. 
The  forward  treatment  of  the  advection  terms  by  ECOM-si  is  quite  dispersive.  The  forward  treatment 
of  the  Coriolis  term  is  unstable  in  that  it  tends  to  cause  the  growth  of  inertial  oscillations.  This  error 
is  small  and  generally  not  noticeable  with  the  small  timesteps  typically  used  in  high-resolution 
coastal  modeling,  but  can  become  significant  when  the  timestep  exceeds  about  200  s.  (The  Coriolis 
term  was  converted  to  an  Adams-Bashforth  treatment  to  avoid  this  timestep  limitation.)  The  for¬ 
ward  time-differencing  scheme  in  ECOM-si  can  also  cause  a  modification  of  the  ambient  stratification 
during  internal  wave  propagation.  This  is  caused  by  numerical  diffusion  due  to  a  phase  (timing) 
error  between  the  vertical  velocity  and  the  temperature  and  salinity  values  in  the  vertical  advection 
terms  of  the  temperature  and  salinity  conservation  equations. 

The  implicit  treatment  of  the  free  surface  in  ECOM-si  and  SZM  is  much  less  accurate  for  the 
propagation  of  surface  waves  than  the  split-explicit  scheme  used  by  POM  in  terms  of  phase  speed 
error  and  damping.  The  phase  speed  errors  for  ECOM-si  and  SZM  are  similar,  but  ECOM-si’s  fully 
implicit  treatment  of  surface  waves  is  significantly  more  damping  than  SZM’s  partially  implicit 
treatment.  We  note,  however,  that  with  the  high  temporal  and  spatial  grid  resolution  typically  used 
in  coastal  modeling,  ECOM-si  and  SZM  generally  simulate  the  tides  fairly  accurately  since  the  long 
wavelength  and  period  of  the  tides  are  well  resolved.  If  coarse  grids  and  large  timesteps  are  used, 
the  accuracy  of  tidal  prediction  with  the  implicit  treatment  of  the  free  surface  may  be  a  concern. 

Sigma  vertical  coordinates  can  suffer  from  problems  with  their  horizontal  advection,  diffusion, 
and  pressure  gradient  terms  in  regions  of  steep  bathymetry.  One  problem  that  can  occur  is  overshoot 
of  the  spatially  centered  advection  term  in  the  bottom  layers  at  a  steep  bottom  slope.  This  is  due 
to  the  sharp  change  (i.e.,  “front”)  that  can  occur  in  the  bottom  sigma  layers  when  a  shallow  point 
is  right  next  to  a  deep  point.  Advection  between  the  shallow  point  and  the  deep  point  can  result  in 
severe  advective  overshoot  due  to  the  large  change  in  the  advected  field  between  the  two  points. 

The  practice  of  subtracting  a  spatially  averaged  profile  from  the  temperature  and  salinity  fields 
when  calculating  horizontal  diffusion  in  sigma  coordinates  can  result  in  significant  spurious  diffusion 
if  the  local  temperature  or  salinity  structure  is  much  different  from  the  mean  profile  that  is  subtracted. 
In  a  downwelling  problem,  the  positive  temperature  anomaly  in  the  downwelling  region,  relative  to 
the  mean  temperature  profile  in  the  model  domain,  resulted  in  a  spurious  warming  of  the  water 
shoreward  of  the  downwelling  area  due  to  diffusion  of  heat  from  the  downwelling  area  toward 
the  shore. 

Z-level  grids  also  have  particular  shortcomings.  If  the  bathymetry  is  truncated  to  the  nearest 
z-level,  as  is  done  with  a  number  of  z-level  models  besides  the  SZM  model  that  is  included  in  this 
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study,  the  problem  that  the  model  actually  solves  is  that  for  the  stair-step  bathymetry  being  used 
in  the  model,  and  not  the  problem  for  the  true  bathymetry  that  the  stair-step  bathymetry  is  approxi¬ 
mating.  (This  may  seem  obvious,  but  there  is  a  tendency  to  think  in  terms  of  the  bathymetry  one 
is  modeling,  rather  than  in  terms  of  the  bathymetry  that  is  actually  in  the  model.) 

On-shore  and  offshore  barotropic  flows  can  be  noticeably  distorted  by  the  stair-step  approximation 
of  a  z-level  grid.  The  horizontal  convergence  of  the  flow  is  focused  at  the  faces  of  the  stair-steps 
rather  than  being  spread  out  over  the  region  of  decreasing  depth,  i.e.,  the  flow  in  the  model  is  the 
flow  that  would  result  if  the  step  were  actually  present.  In  a  problem  with  an  on-shore,  barotropic 
tidal  flow,  the  isotherms  were  distorted  by  the  vertical  “jets”  that  occurred  at  the  faces  of  the  steps. 
In  the  real  ocean,  internal  waves  are  generated  by  such  tidal  flows  in  regions  where  there  are  sharp 
bathymetric  changes,  but  this  should  not  happen  in  a  region  where  the  bathymetry  is  changing 
gradually.  Topographic  shelf  waves,  which  depend  on  changes  in  bottom  depth  for  their  existence, 
can  be  quite  distorted  by  a  stair-step  approximation  to  a  smoothly  varying  bathymetry.  The  best 
solution  to  these  problems  with  z-level  vertical  coordinates  may  be  to  truncate  the  bottom  grid  cell 
of  the  z-level  grid  to  the  bathymetry,  rather  than  rounding  the  bottom  depth  to  the  nearest  model 
level.  This  adds  complication,  but  is  being  done  in  some  z-level  models. 

The  Crank-Nicolson  scheme  used  for  vertical  mixing  in  SCRUM  2.1,  where  the  values  being 
diffused  are  evenly  weighted  at  the  old  and  new  time  levels,  causes  diffusive  overshooting  if  the 
vertical  diffusion  coefficients  are  large,  which  they  frequently  are  in  realistic  mixing  situations.  The 
result  is  that  vertical  gradients  in  the  field  being  diffused  (temperature,  salinity,  or  velocity)  can  be 
reversed  (or  distorted)  between  adjacent  vertical  points.  The  computation  does  not  blow  up,  but  the 
resulting  mixing  can  be  noisy  and  inaccurate.  All  the  other  models  use  fully  implicit  vertical 
mixing,  with  the  field  being  diffused  weighted  fully  at  the  new  time  level,  to  avoid  this  problem. 

Checkerboard  mixing,  where  a  fluctuation  in  the  mixed-layer  depth  sets  up  at  alternate  gridpoints, 
was  found  to  occur  with  the  models  under  certain  conditions  of  light  to  moderate  winds  and  surface 
heating  (or  a  positive  surface  buoyancy  flux).  The  checkerboard  mixing  occurs  because  of  the 
horizontal  averaging  that  is  used  when  computing  vertical  mixing  on  a  C  grid,  where  the  velocity 
and  the  temperature-salinity  points  are  at  different  locations.  As  a  practical  matter,  checkerboard 
mixing  is  not  usually  seen  in  realistic  simulations  because  of  the  temporal  and  spatial  changes  in 
the  surface  forcing  that  tend  to  suppress  or  mask  the  phenomenon. 
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Appendix  A 

STABILITY  ANALYSIS  FOR  SURFACE  WAVE  PROPAGATION 


The  properties  of  a  numerical  scheme  can  be  determined  by  performing  a  stability  analysis  of 
the  scheme.  Here  we  will  analyze  the  propagation  of  surface  waves  in  the  models  using  a  procedure 
referred  to  as  a  von  Neumann  stability  analysis.  This  analysis  can  provide  information  about  the 
stability  and  damping  properties  of  the  scheme  and  about  the  dispersion  or  phase  speed  error.  A 
good  discussion  of  the  von  Neumann  stability  analysis  can  be  found  in  Mesinger  and  Arakawa 
(1976),  and  a  discussion  of  the  analysis  of  the  Asselin  filter  can  be  found  in  Asselin  (1972).  We 
will  look  at  the  propagation  of  surface  waves  for  the  case  where /=  0,  which  significantly  simplifies 
the  analysis  of  the  numerical  schemes. 

The  general  procedure  is  to  substitute  a  solution  of  the  form 

Ux,  t)=Z(t)eikx  (Al) 

into  the  linearized  numerical  equations,  where  Z  is  the  time-dependent  part  of  the  solution  and  k 
is  the  horizontal  wavenumber.  This  allows  the  determination  of  the  behavior  of  a  single  Fourier 
spatial  mode  of  wavelength  L  =  2n/k.  One  then  solves  for  the  amplification  factor  X  =  Z(n  +  lVz^n\ 
which  is,  in  general,  complex.  The  magnitude  or  modulus  of  X  denotes  the  stability  and  damping 
properties  of  the  scheme,  i.e.,  how  the  magnitude  of  the  solution  changes  on  successive  timesteps. 
From  the  phase  change  of  X  each  timestep,  the  propagation  speed  of  the  waves  in  the  numerical 
solution  can  be  computed. 

The  linearized  equations  for  the  propagation  of  surface  gravity  waves  in  the  x-direction  for 
/  =  0  are 


du 

(A2) 

Jt  = 

dx 

dt 

du 

(A3) 

~dt~ 

~h12- 

Al.  STABILITY  ANALYSIS  FOR  SURFACE  WAVES  IN  POM 

For  POM,  the  numerical  scheme  for  the  propagation  of  surface  waves  is  an  explicit  leapfrog 
scheme  with  an  Asselin  temporal  filter.  The  finite  difference  forms  of  (A2)  and  (A3)  for  this 
numerical  scheme  are 


.  i  ~Q.  i 


u(n  +  1)  =  A[u(n~1)]-2Atg- 


l  +  r 


l~2 


Ax 


(A4) 
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(n)  («) 

«.  1  -“._1 

^(«  + 1)  =  A  [^(^  - 1)]  _  2A  tH  <  +  2a^  *  2  .  (A5) 

A  [  ]  denotes  application  of  the  Asselin  filter,  which  is  defined  as 

A [«(”)]  =  +  v(m^w  ~ -  2k^  +  +  *))  ,  (A6) 

where  v  is  the  filter  coefficient. 

Substitution  of  the  solution  form  (Al)  into  the  numerical  equations  (A4)-(A6)  and  solving  for 
X  gives 

\  =  (iq  +  v)±((l-v)2-q2)2,  (A7) 

where 

.  sl2Ar  .  kAx „  ,  A  „ 

<l  =  (gH) 2~^  sin  (— )  .  (A8) 

Hence,  the  amplification  each  timestep  is 

|X|  =  ((1  -  v)2  +  v2  +  2v  ((1  -  v)2  -  q2)*  )2  (A9) 

and  the  phase  change  each  timestep  is 

coAr  =  tan_1  (<?/((  1  -  v)2  -  q2^2  +  v) ,  (A10) 

where  co  is  the  frequency.  The  phase  speed  is  then  given  by 

C  =  j  =  ^tan_  *($  /  (((!  -  V)2  -  ^ 2)2  +  v))  .  (Al  1) 

Note  that  for  v  =  0  we  get  |a|  =  1,  i.e.,  there  is  no  damping  of  the  basic  leapfrog  scheme  without 
the  Asselin  filter.  The  damping  is  due  entirely  to  the  filter. 


A2.  STABILITY  ANALYSIS  FOR  SURFACE  WAVES  IN  ECOM-SI 

For  ECOM-si,  the  numerical  scheme  for  the  propagation  of  surface  waves  is  a  two-time-level, 
fully  implicit  scheme,  i.e.. 


>  +  l)  ,.(»  +  !) 

<=>.  .  l  -<=>.  l 


I  +: 


•A  tg- 


l~2 


Ax 


(A12) 


A  Comparison  of  Several  Coastal  Ocean  Models 


87 


( n  +  1)  (n  +  1) 

u  l  -«._i 

+ 1)  _  ^(n)  _  A  tH - !_ - —  .  (A1 3) 

Again,  substitution  of  the  solution  form  (Al)  into  the  numerical  equations  (A12)-(A13)  and  solv¬ 
ing  for  X  gives 


X  = 


1  ±  iq 


1  +  9 


2  ' 


(A14) 


Hence,  the  amplification  each  timestep  is 

N  =  (l  +q2)~1i  (A15) 

and  the  phase  change  each  timestep  is 

o)A t  =  tan-1^) .  (A16) 

The  phase  speed  is  then  given  by 

<A17> 

Note  that  the  amplification  factor  is  always  less  than  1,  i.e.,  the  scheme  is  always  stable. 


A3  STABILITY  ANALYSIS  FOR  SURFACE  WAVES  FOR  SZM 

For  SZM,  the  numerical  scheme  used  in  this  report  for  the  propagation  of  surface  waves  is  a 
trapezoidal  implicit,  leapfrog  scheme  with  an  Asselin  temporal  filter.  The  finite  difference  form  of 
(A4)-(A5)  for  this  scheme  is 


An  + 1)  (n  +  1) 
Q.  l  -Q.  l 


u(n  +  V=A[u(n~V]-Atg- 


(n  -  1)  (n  - 1) 

.  .  .  m  i  ]-a[£  i  j 

1  +  2  ‘-2  A#  l+2  l~2 
-A  tg - 


Ax 


Ax 


(A18) 


(n  +  1)  (n  +  1)  (n-1)  (n-1) 

U.  1  -M.  1  A[ M  l  ]-A[m._i  j 

t(”  +  1)=A[^"~1>]-ArJj'  —  .  *  - - A  tH-  1  2  2 


Ax 


Ax 


(A19) 


Substitution  of  the  solution  form  (Al)  into  the  numerical  equations  (A18),  (A19),  and  (A6) 
(for  the  Asselin  filter)  and  solving  for  X  gives 


X  =  (v  ±  ((1  -  v)2  +•  q2  (1  -  2v))2)  -  +  I^f  . 

1+q2 


(A20) 
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Hence,  the  amplification  each  timestep  is 


|X|  =  (v  +  ((1  -  v)2  +  q2  (1  -  2v))2)(1  +  q2)  2 
and  the  phase  change  each  timestep  is 


(A21) 


toA  t  =  tan  ^(q) ,  (A22) 

where  od  is  the  frequency.  The  phase  speed  is  then  given  by 

c  =  I  =  A7Tlan‘,(«)'  (A23) 

The  amplification  is  always  less  than  or  equal  to  1  and  the  scheme  is  always  stable.  When  v  =  0, 
we  get  |X|  =  1,  i.e.,  there  is  no  damping  of  the  waves.  Hence,  like  POM,  the  damping  of  the  waves 
is  due  entirely  to  the  Asselin  filter.  The  phase  speed  error  is  the  same  as  for  ECOM-si,  and  is  not 
affected  by  the  value  of  v. 

Figures  37-39  show  plots  of  the  phase  speed  error  and  damping  of  the  numerical  schemes  as 
a  function  of  the  Courant  Number  Cu  =  (gH)''2  At/ Ax  and  the  number  of  gridpoints  per  wavelength 
that  are  used  to  resolve  the  wave.  The  phase  speed  error  is  expressed  as  the  percent  error  relative 
to  the  analytically  determined  phase  speed  for  the  wave,  and  the  damping  is  expressed  as  the 
percent  damping  per  wave  period  based  on  the  analytically  calculated  wave  period.  The  damping 


Fig.  37  —  (a)  Phase  speed  error  and  (b)  damping  for  the  numerical  scheme  used  by  POM  for  the  free  surface  mode. 
The  phase  speed  error  is  expressed  as  the  percent  error  relative  to  the  analytically  calculated  phase  speed.  The  damping 
is  expressed  as  the  percent  damping  per  wave  period  based  on  the  analytically  calculated  wave  period.  The  Asselin 
filter  coefficient  for  this  calculation  is  set  to  0.05.  Contour  labels  less  than  0.1  are  multiplied  by  10,000. 
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(a)  (b) 


GRIDPOINTS  PER  WAVELENGTH 


Fig.  38  —  Same  as  Fig.  37,  but  for  ECOM-si 


error  for  POM  and  SZM  depends  on  the  value  of  the  Asselin  filter  coefficient  v.  For  the  plots  in 
Figs.  37  and  39,  v  was  set  to  the  typically  used  value  of  0.05.  Values  are  shown  for  POM  only 
for  values  of  Cu  less  than  0.5,  since  POM’s  explicit  scheme  is  unstable  for  Cu  >  0.5. 

The  phase  speed  error  for  all  three  schemes  is  the  same  at  very  low  values  of  Cu,  but  whereas 
the  phase  speed  error  for  POM  decreases  as  Cu  increases  (up  to  the  point  at  which  POM’s  explicit 
scheme  becomes  unstable  at  Cu  =  0.5),  the  phase  speed  error  for  ECOM-si  and  SZM  (which  is  the 
same)  increases.  The  damping  for  POM  and  SZM  is  similar  for  Cu  <  0.5  and  is  much  less  than 


COURANT  NUMBER 
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the  damping  of  ECOM-si  at  the  same  value  of  Cu.  Of  course,  it  must  be  remembered  that 
ECOM-si  and  SZM  will  typically  be  using  much  larger  values  of  Cu  for  surface  waves  than  POM 
because  of  the  small,  barotropic  timestep  used  by  POM. 


(a)  (b) 


GRIDPOINTS  PER  WAVELENGTH 


Fig.  39  —  Same  as  Fig.  37,  but  for  SZM.  The  Asselin  filter  coefficient  for  this  calculation  is  set  to  0.05. 


Appendix  B 

INTERNAL  WAVE  PROPAGATION  ERROR  WITH  ECOM-SI 


The  propagation  of  internal  gravity  waves  with  ECOM-si  resulted  in  some  modification  of  the 
ambient  stratification.  An  analysis  of  internal  gravity  wave  propagation  with  ECOM-si ’s.  forward 
time-differencing  scheme  led  to  the  conclusion  that  this  is  primarily  due  to  numerical  diffusion  of 
the  ambient  stratification  caused  by  a  phase  (timing)  error  between  the  salinity  and  velocity  terms 
in  the  salinity  equation’s  vertical  advection  term. 

The  main  terms  governing  the  propagation  of  internal  waves  are  the  horizontal  baroclinic 
pressure  gradient  and  the  vertical  advection  of  density.  For  internal  wave  propagation  as  described 
in  Sec.  6.0,  i.e.,  for  propagation  in  the  ^-direction  with  the  vertical  density  gradient  being  due  to 
salinity  stratification,  and  with  /=  0,  the  relevant  momentum  and  salinity  equations  are 


du  _  1  dp 

dt  po  dx  ’ 


(Bl) 


dS  d(wS ) 

dt  dz 


(B2) 


The  salinity  field  can  be  split  into  two  parts,  one  part  S0  that  describes  the  ambient  stratification, 
which,  for  the  problem  being  considered  here,  depends  only  on  depth,  and  another  part  S'  that 
describes  the  changes  in  the  ambient  salinity  field  due  to  the  propagating  wave,  i.e., 

S(x,  z,  t)  =  S0(z)  +  S'  (x,  z,  t). 

Hence,  the  salinity  equation  can  be  written  as 

dS  HwSq)  d(wS') 

dt  dz  dz 


(B3) 


(B4) 


For  small  amplitude  waves,  the  main  contribution  to  the  vertical  advection  of  salinity  is  due  to  the 
vertical  advection  of  the  ambient  stratification  SQ.  The  contribution  from  the  vertical  advection  of 
S'  is  small. 
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With  the  two-time-level  numerical  scheme  used  by  ECOM-si,  the  momentum  is  first  updated 
using  the  horizontal  pressure  gradient  that  has  been  calculated  using  the  old  temperature  and  salinity 
values 


?o  dx 


(«) 


(B5) 


and  then  the  new  salinity  is  calculated  using  the  newly  calculated  velocities  for  the  advection  term 

(B6) 


dz 


This  calculation  sequence  can  be  conceptualized  as  a  leapfrogging  of  the  momentum  and  density 
terms  in  time  (Fig.  40),  even  though  the  model  is  not  explicitly  written  as  if  this  were  the  case.  With 
this  leapfrogging,  the  velocity  and  salinity  fields  are  offset  from  each  other  by  half  a  timestep  and 
the  main  terms  describing  the  internal  wave  propagation  are  accurately  centered  in  time,  i.e.,  the 
salinity  used  for  the  calculation  of  the  horizontal  pressure  gradient  is  located  between  the  old  and 
new  values  of  u,  and  the  velocity  used  for  advection  of  salinity  is  located  between  the  old  and  new 
values  of  5. 

A  quantity,  however,  that  is  not  centered  in  time  is  the  value  of  S  used  in  the  advection 
term.  The  value  that  is  used  is  the  old  value  S(n\  which  effectively  lags  w^n  +  U  by  half  a  timestep. 
Since  most  of  the  contribution  to  the  salinity  advection  is  from  the  ambient  salinity  stratification 
Sa,  only  a  secondary  contribution  to  the  salinity  advection  suffers  from  a  phase  lag  error.  However, 
it  is  this  phase  lag  of  S  in  the  advection  term  that  results  in  the  modification  of  the  ambient  salinity 
stratification. 

If  (B4)  is  horizontally  averaged  over  one  wavelength,  we  get  an  expression  for  the  rate  of 
change  with  time  of  the  horizontally  averaged  salinity  S(z,t) 


dS  d(wS0)  d(wS') 

dt  dz  dz 


(B7) 


where  the  overbar  indicates  a  horizontal  average.  The  first  term  on  the  right  side  of  (B7)  is  zero,  since 
S0  has  no  horizontal  variation  and  w  =  0.  The  second  term  on  the  right  side  of  (B7)  describes  the  time 


UM  v(n>  w(n)  u{n+1)  v<rt+1)  w(/l+1) 
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Fig.  40  —  Schematic  showing  the  effective  leapfrogging 
in  time  of  the  momentum  and  salinity  conservation 
equations  of  ECOM-si  for  the  main  internal  gravity  wave 
terms 
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rate  of  change  of  the  horizontally  averaged  salinity  due  to  the  vertical  divergence  of  the  salt 
flux  defined  by  wS'.  Using  (74)  and  (75)  from  Sec.  6.0,  we  can  write  vv5'  as 


yki  2CO  o 

wS’  =A  — —  sin  ( mz )  sin  (kx  -  tut)  cos  (kx-tot)  . 


(B8) 


If  (B8)  is  horizontally  averaged  over  one  wavelength,  the  product  of  sin  (kx  -  w t)  and  cos  (kx  -  w t) 
integrates  to  zero,  i.e.,  there  is  no  net  vertical  flux  of  salt  due  to  the  passage  of  the  internal  wave. 
However,  if  there  is  an  error  in  the  phase  of  S'  relative  to  w,  the  horizontal  average  of  (B8)  will 
not  be  zero. 


Since  we  estimate  that  for  internal  wave  propagation  in  ECOM-si,  the  value  of  S  in  the  salinity 
advection  term  effectively  lags  w  by  half  a  timestep,  we  get 


oN  oo  2 

h>5'  =  A  sin  (mz)  sin  (kx -  cat)  cos  (kx-  <x>(t- A  t/2))  . 


(B9) 


Using  some  trigonometric  relations,  this  becomes 


2 

*y j^J~  Q)  ry  ry  £ 

wS'  =  -A  sin  (mz)  sin  (kx -  tot)  sin  (o>y)  . 


Taking  the  horizontal  average  over  one  wavelength  gives 


(BIO) 


—  oiV  oo  9  A  t 

wS  =  -A Z~Xp—  sin  (mz)  sin  (to—) .  (B 1 1) 

Zpg  2 

Hence,  the  mean  vertical  flux  of  salinity  due  to  the  phase  lag  of  S'  with  respect  to  w  is  downward. 
Note  that  the  direction  of  this  flux  is  against  the  ambient  salinity  gradient.  The  rate  of  change  of 
the  salinity  with  time  is  given  by  the  negative  of  the  spatial  derivative  of  the  flux  and  is 


dS 

dt 


d(wS') 

dz 


oA j2u>m  ...  A  t. 

=A  — — sm  (2m  z)  sin  (<o— ) 
2P g  2 


(B12) 


The  results  of  the  analysis  predict  that  the  rate  of  change  of  the  background  stratification  with  time 
depends  on  the  sine  of  the  phase  error  between  w  and  S  in  the  salinity  advection  term,  and  is 
proportional  to  the  square  of  the  amplitude  of  the  internal  wave. 

Figure  41a  shows  a  comparison  of  5'  calculated  by  (B12)  with  the  result  from  ECOM-si  at  24  h 
for  the  propagation  of  the  8-km  wavelength  internal  wave  reported  in  Sec.  6.0  with  At  =  200  s.  Also 
shown  are  comparisons  from  two  similar  experiments  -  one  in  which  the  amplitude  of  the  wave  was 
increased  from  100  to  200  cm  with  At  kept  at  200  s  and  another  in  which  the  amplitude  was  increased 
to  400  cm  and  the  timestep  was  increased  to  400  s.  The  comparisons  show  good  agreement  between 
the  analysis  and  the  results  from  ECOM-si  in  terms  of  the  structure  and  amplitude  of  the  change 
in  the  ambient  stratification.  At  the  higher  wave  amplitudes,  the  analysis  appears  to  overestimate 
the  change  in  the  ambient  stratification  slightly. 
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CHANGE  (psu) 


Fig.  41  —  Change  in  ambient  stratification  during  internal  wave  propagation  for  ECOM-si 
(solid  line)  and  as  predicted  by  numerical  analysis  (dashed  line)  after  24  h.  The 
wavelength  is  8  km.  The  internal  wave  amplitude  and  timestep  are  (a)  100  cm  and 
200  s,  (b)  200  cm  and  200  s,  and  (c)  400  cm  and  400  s. 


Appendix  C 


ANALYTICAL  SOLUTION  FOR  TOPOGRAPHIC  WAVES 
WITH  STAIRSTEP  BOTTOM 


The  linearized  equations  for  barotropic  shelf  waves  propagating  in  the  x-direction  along  a 
straight  coast  with  bathymetry  varying  only  in  the  offshore  direction  are 


du  dt, 

(Cl) 

dv  dt 

(C2) 

dx  dy 

(C3) 

Note  that  the  time  derivative  has  been  scaled  out  of  the  continuity  equation.  This  term  is  not 
essential  to  the  physics  of  small-amplitude,  topographic  waves,  and  its  elimination  simplifies  the 
solution  of  the  equations  (LeBlond  and  Mysak  1978). 

Cross-differentiation  of  the  momentum  equations  yields  the  vorticity  equation 


d  dv  d«  du  dv 

dt  dx  dy  dx  +  dy  ' 


(C4) 


Based  on  the  continuity  equation,  we  can  define  a  stream  function  W  where 


1  dW 

U~~H  dy  ’ 


(C5) 


1  dW 
V  “  Hlx 


(C6) 


Substituting  the  stream  function  into  the  vorticity  equation,  we  get 


1  d2W 


±  v  ,  d_  I3f  _/_aff3V 
H  qx2  +  dy  W  dy  2  dy  dx 


(C7) 
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We  can  assume  a  traveling  wave  solution  in  x  of  the  form 


(C8) 


where  k  is  the  along-shore  wavenumber  and  co  is  the  frequency. 

Substituting  into  (C7)  and  dropping  the  prime  on  W,  we  get  an  equation  for  the  cross-shore 
structure  of  the  topographic  wave 


d  1  dW  fk  dH 
dy{H  dy)  +  \,H2  dy 


H 


)  'P  =  0  . 


This  equation  can  be  solved  for  a  particular  bathymetry  H(y). 
For  a  region  of  constant  depth,  (C9)  simplifies  to 


(C9) 


^4-)fc2W  =  0.  (CIO) 

dy2 


Solutions  of  this  equation  have  the  form  of  cosh  and  sinh  functions. 

Since  a  step-wise  bathymetry  consists  of  a  series  of  adjoining  regions  of  uniform  depth,  with 

a  discontinuity  in  the  depth  between  each  region,  the  solution  can  be  represented  as  cosh  or  sinh 

functions  within  each  region  of  uniform  depth,  with  the  amplitude  of  the  functions  determined  by 
appropriate  matching  conditions  between  the  regions. 

For  a  domain  extending  from  y  =  0  to  y  =Ly,  in  which  the  depth  between  y;-_  1  and  yy  is  Hj  for 
j  =  1  to  j  =  m,  the  solution  within  the  region  from  y.-  _  i  to  yy  can  be  represented  as 

Wj  =  Aj  sinh(Ay)  +  Bj  cosh(Ay) .  (Cll) 

Boundary  conditions  of  zero  transport  (W  =  0)  at  y  =  0  and  y  =  Ly  give 

=  sinh(Ay)  ,  (C12) 

=  Am  sinh(£(y  -  Ly )) .  (C13) 

The  unknowns  Ay,  Bp  and  to  can  be  determined  from  the  matching  conditions  between  the 
regions  at  each  yj  (Stocker  and  Hutter  1987) 

Wy  =  Wy  +  1,  (C14) 
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H,  v  dy 


co  . J  +  V 


Hj  + 1 


(C15) 


This  coupled  set  of  2m -2  equations  for  the  cross-shore  structure  of  the  topographic  wave  for  a 
stairstep  bottom  can  be  solved  numerically  using  a  shooting  method. 


